0

三重対角行列を使用した 3 次スプライン補間について、MATLAB での宿題のためにこのコードを書きました。私はアルゴリズムのすべてのステップに従いますが、実際にはエラーが見つかりません。コードは 2 年生の関数では問題ありませんが、たとえば、sin(x) を配置すると、結果はスプラインではなく、他の関数では問題がないため理由がわかりません。誰でもエラーを見つけるのを手伝ってもらえますか? ありがとう

3 次スプライン スクリプト:

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

三重対角行列:

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

tk ポイントでのスプライン評価:

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;
4

1 に答える 1

0

それは、Matlab をfor正しく使用していないためです。

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

最後のfor読み取りfor i = n - 1: 1が実行されない場合は、次のように記述します。

for i = n - 1:-1:1

その後、動作します。sin(x) だけでなく、以前の試みでも機能しなかったことに注意してください。

ここに画像の説明を入力

于 2014-11-07T11:53:42.010 に答える