Сплайн кубический с тридиагональной матрицей

Question

Я написал этот код для своей домашней работы на MATLAB о кубической сплайновой интерполяции с трехдиагональной матрицей. Я следую всем шагам алгоритма, и я действительно не нахожу свою ошибку. Код в порядке со второстепенными функциями, но когда я ставлю, например, sin (x), результат не является сплайном и не знает, почему, потому что с другой функцией у меня нет проблем. Может ли кто-нибудь помочь мне найти ошибку? Благодаря

кубического сплайна SCRIPT:

close all 
clear all 
clf reset 

ff = @(x) sin(x); 
x = [-2 0 2 4 6 7 8 9 10 11 12]; 
for i = 1: length(x), 
    omega(i) = ff(x(i)); 
end 

n = length(x); 
h = zeros(n - 1, 1); 
for i = 1: n - 1, 
    h(i) = x(i + 1) - x(i); 
end 

a = zeros(n, 1); 
b = zeros(n, 1); 
d = zeros(n, 1); 
f = zeros(n, 1); 

for j = 2: n - 1, 
    a(j) = 2*(h(j) + h(j - 1)); 
    b(j) = h(j - 1); 
    d(j) = h(j); 
    f(j) = 6 * (omega(j + 1) - omega(j)/h(j)) - (6 * (omega(j) - omega(j - 1))/h(j - 1)); 
end 

% Starting conditions 
a(1) = -2; 
f(1) = 0; 
a(n) = 4; 
f(n) = 0; 

% Coefficents 
c = tridiag(a, b, d, f); 

t = linspace (x(1), x(n), 301); 
for k = 1: length(t), 
    tk = t(k); 
    y(k) = spline_aux(x, omega, c, tk); 
    z(k) = ff(tk); 
end 

plot(t, z, 'b-', 'linewidth', 2) 
hold on; 
plot(t, y, 'r.', x, omega, 'go') 
grid on 

TRIADIAGONAL MATRIX:

function [x] = tridiag(a, b, d, f) 

n = length(f); 
alfa = zeros(n, 1); 
beta = zeros(n, 1); 

alfa(1) = a(1); 
for i = 2: n, 
    beta(i) = b(i)/alfa(i - 1); 
    fprintf (' i: %d  beta: %12.8f\n', i, beta(i)) 
    alfa(i) = a(i) - (beta(i)*d(i - 1)); 
    fprintf (' i: %d  alfa: %12.8f\n', i, alfa(i)) 
end 

y(1) = f(1); 
for i = 2: n, 
    y(i) = f(i) - beta(i)*y(i - 1); 
end 

x(n) = y(n)/alfa(n); 
for i = n - 1: 1, 
    x(i) = (y(i) - (d(i)*x(i + 1)))/alfa(i); 
end 

СПЛАЙН ОЦЕНКА В Т.К. POINT:

function [s] = spline_aux(x, w, c, tk) 

n = length(x); 

h = zeros(n - 1, 1); 
for i = 1: n - 1, 
    h(i) = x(i+1) - x(i); 
end 

for i = 1: n - 1, 
    if (x(i) <= tk && tk <= x(i+1)) 
     break 
    end 
end 

s1 = c(i)*((x(i+1) - tk)^3)/(6*h(i)); 
s2 = c(i+1)*((tk - x(i))^3)/(6*h(i)); 
s3 = (w(i)/h(i) - (c(i)*h(i)/6))*(x(i+1) - tk); 
s4 = (w(i+1)/h(i) - (c(i+1)*h(i)/6))*(tk - x(i)); 
s = s1 + s2 + s3 + s4; 

Answer 1

Это потому, что вы не используете Matlab-х for правильно

в функции function [x] = tridiag(a, b, d, f)

последний for читает for i = n - 1: 1 но никогда не будет выполнять, вы shoudl writte:

for i = n - 1:-1:1

Затем работает. Следует заметить, что didt работу в любой из предыдущих попыток, не только грех (х)