chap1_3.m

利用MATLAB解决<偏微分方程数值算法>中的一些问题,如插值求解

function[Y,E,e]=CHAP1_3(m,Z)%m为等分的分数,Z为欲知该点函数值的点集,Y为对应的函数值.E为对应的误差绝对值,e为最大误差．
for i=1:m-1
    X(1,i)=i*pi/(2*m);
    B(1,i)=2+(pi/(2*m))^2*(X(1,i)-0.5)^2;
    D(1,i)=(pi/(2*m))^2*(X(1,i)^2-X(1,i)+1.25)*sin(X(1,i));
end
for i=1:m-2
    A(1,i)=-1;
end
for i=1:m-2
    R(1,i)=-1;
end
D(1,1)=D(1,1)+0;
D(1,m-1)=D(1,m-1)+1;
U=Thomas2(A,B,R,D);%以上部分求出了截点处函数值．
P=U';%由于求出的U为一个列矩阵;故转置．
Q(1,1)=0;W(1,1)=0;
for i=2:m
    Q(1,i)=P(1,i-1);W(1,i)=X(1,i-1);
end
Q(1,m+1)=1;W(1,m+1)=pi/2;%Ｗ为包括０，派的截点，Ｑ为对应函数值．
for i=1:length(Z)
    k=1;
while(Z(1,i)>W(1,k))
k=k+1;
end
Y(1,i)=Q(1,k-1)+((Q(1,k)-Q(1,k-1))/(W(1,k)-W(1,k-1)))*(Z(1,i)-W(1,k-1));
end%利用线形函数性质求给定点处的函数值．
for i=1:length(Z)
    E(1,i)=abs(Y(1,i)-sin(Z(1,i)));
end
for i=1:m+1
    E1(1,i)=abs(Q(1,i)-sin(W(1,i)));
end
e=max(E1);
plot(W,Q);hold on;
plot(W,E1,'r');
q=0:0.01:pi/2;
p=sin(q);
plot(q,p,'r-');