suanli1_1_3.m

孙志忠编写的《偏微分方程数值解法》一书中例题1.1.3的程序。用紧差分格式解常微分方程边值问题。

%算例1.1.3，见课本第17页
%首先网格剖分
xl=0.;
xr=pi;%左右端点
n=20;
h=(xr-xl)/n;%步长
%解差分格式形成的方程组
U=zeros(n-1,1);
aa=-1+h^2/12;
bb=2+5/6*h^2;
A=diag(bb*ones(n-1,1))+diag(aa*ones(n-2,1),1)+diag(aa*ones(n-2,1),-1);%系数矩阵
F=zeros(n-1,1);
for i=1:n-1
    F(i)=h^2/12*(exp(xl+(i-1)*h)*(sin(xl+(i-1)*h)-2*cos(xl+(i-1)*h))+10*exp(xl+i*h)*(sin(xl+i*h)-2*cos(xl+i*h))+exp(xl+(i+1)*h)*(sin(xl+(i+1)*h)-2*cos(xl+(i+1)*h)));%右端项
end
U=inv(A)*F;%得到数值解
%精确解与数值解误差绝对值的最大值
ana=zeros(n-1,1);
for i=1:n-1
    ana(i)=exp(xl+i*h)*sin(xl+i*h);%精确解
end
wucha=max(abs(U-ana))
%可视化处理
x=xl:h:xr;
U=[0 U' 0];%U表示数值解
ana=[0 ana' 0];%ana表示精确解
plot(x,U,'r*',x,ana);