cubic_spline.m

来自「该程序包含矩阵最大特征值上、下界估计、线性方程组的求解、非线性方程组的一个最新迭」· M 代码 · 共 43 行

%本程序是用来给出奇异扰动边值问题的三次样条方法，结合三次样条函数后的方法，精度明显提高了血多！


clear
format long
e=1/10^(0);
h=1/32;
nh=1/h;
p=1;
a=p*e-h^2*p^2/12;
b=h^2*p/12-e;
c=e*h^2/12;
d=e^2;
for i=1:nh+1
    f(i)=-cos(pi*(i-1)*h)*cos(pi*(i-1)*h)-2*e*pi^2*cos(2*pi*(i-1)*h);
    f1(i)=-2*sin(pi*(i-1)*h)^2*pi^2+2*cos(pi*(i-1)*h)^2*pi^2+8*e*pi^4*cos(2*pi*(i-1)*h);
    f2(i)=-8*cos(pi*(i-1)*h)^2*pi^4+8*sin(pi*(i-1)*h)^2*pi^4-32*e*pi^6*cos(2*pi*(i-1)*h);
end
y(1)=0;
y(nh+1)=0;
for j = 1:nh-1
    if  j ~= nh-1
        A(j,j) =-12*d/h^2-4*a; ;
        B(j,j)=4;
        A(j+1,j) =6*d/h^2-a;
       B(j+1,j)=1;
        A(j,j+1) =6*d/h^2-a;
       B(j,j+1)=1;
    else
        A(j,j) =-12*d/h^2-4*a;
      B(j,j)=4;
    end
end
    D(1:nh-1)=b*B*f(2:nh)'+c*B*f1(2:nh)';
    D(1)=D(1)+b*f(1)-(6*d/h^2-a)*y(1)+c*f1(1);
    D(nh-1)=D(nh-1)-(6*d/h^2-a)*y(nh+1)+b*f(nh+1)+c*f1(nh+1);
    y(2:nh)= lufact(A,D');
for i=1:nh+1
    u(i)=(exp(-(1-(i-1)*h)/e^(1/2))+exp(-(i-1)*h/e^(1/2)))/(1+exp(-1/e^(1/2)))-cos(pi*(i-1)*h)*cos(pi*(i-1)*h);
end
x=0:h:1;
plot(x,y,'--',x,u,'+')
max(abs(y-u))