sample.m

来自「三次样条插值函数的三弯矩法」· M 代码 · 共 34 行

function [M,S]=sample()
x=linspace(-1,1,11);
m=-18./(1+9.*x.^2).^2.*x;
f1=1./(1+9.*x.^2);
f2=zeros(1,10);f3=zeros(1,9);
%equal inervals,get coeffitientmatrix to get thesolution A
h=0.2000;d=zeros(1,11);
lamda=0.5;miu=0.5;
for i=1:10
    f2(1,i)=(f1(1,i+1)-f1(1,i))./h;
end
for i=1:9
    f3(1,i)=(f2(1,i+1)-f2(1,i))./2.*h;
end
d(1,2:end-1)=6.*f3;
d(1,1)=6./0.2000.*(f2(1,1)-m(1,1));
d(1,end)=6./0.2000.*(m(1,end)-f2(1,end));
A=zeros(11);
A(1,1)=2;A(1,2)=1;
A(11,10)=1;A(11,11)=2;
for i=1:9
    A(i+1,i:i+2)=[miu 2 lamda];
end
%SOR try
M=SOR(A,d',1,ones(11,1));
syms  X;
for i=1:10
    S(i,1)=expand((x(i+1)-X).^3./6.*0.2.*M(1,i)+(X-x(i)).^3./6.*0.2.*M(1,i+1)+(f1(1,i)-0.04.*M(1,i)./6).*(x(1,i+1)-X)./0.2+(f1(1,i+1)-0.04.*M(1,i+1)./6).*(X-x(1,i))./0.2);
    w=(0.2.*(i-6):0.01:0.2.*(i-5));
    y=subs(S(i,1),X,w);
    hold on;
    plot(w,y);
end