# 偏微分方程数值解法 上机实验(3.12)

0.5,0.0 0.5,0.1 0.5,0.2 0.5,0.3 0.5,0.4 0.5,0.5 0.5,0.6 0.5,0.7 0.5,0.8 0.5,0.9 0.5,1.0

0.5,0.0 0.5,0.1 0.5,0.2 0.5,0.3 0.5,0.4 0.5,0.5 0.5,0.6 0.5,0.7 0.5,0.8 0.5,0.9 0.5,1.0

0.5,0.0 0.5,0.1 0.5,0.2 0.5,0.3 0.5,0.4 0.5,0.5 0.5,0.6 0.5,0.7 0.5,0.8 0.5,0.9 0.5,1.0

4.数值解的最大误差

h (1/10) (1/20) (1/40) (1/80) (1/160)

v (1/10) (1/20) (1/40) (1/80) (1/160)

Emax E(2h,2v)/E(h,v) 0.0050206 * 0.0012915 3.88728065 3.27E-04 3.949675841 8.24E-05 3.968446602 2.07E-05 3.980676329

6.用表中数值进行外推，由于是 O（h^2+v^2）的收敛阶 故 u=(4u2-u1)/3

0.5,0.0 0.5,0.1 0.5,0.2 0.5,0.3 0.5,0.4 0.5,0.5 0.5,0.6 0.5,0.7 0.5,0.8 0.5,0.9 0.5,1.0

xn=0:h:1; tn=0:v:1; n=(1/h); m=(1/v); %初始化边界条件 for i=1:n+1 u(1,i)=exp(xn(i))*sin(1/2); end for i=1:m+1 u(i,1)=sin((1/2)-tn(i)); end for i=1:m+1 u(i,n+1)=exp(1)*sin((1/2)-tn(i)); end %构建方程组 for i=1:n-1 for j=1:n-1 if j==(i-1)||j==(i+1) q(i,j)=-r/2; else if i==j q(i,j)=1+r; end end end end for i=1:n-1 for j=1:n-1 if j==(i-1)||j==(i+1) p(i,j)=r/2; else if i==j p(i,j)=1-r; end end end end for k=2:m+1 for j=1:n-1 if j==1 g(j)=(r*(u(k-1,1)+u(k,1))/2+v*f(xn(1+1),(tn(k-1)+tn(k)/2))); else if j==n-1 g(j)=(r*(u(k-1,n+1)+u(k,n+1))/2+v*f(xn(j+1),((tn(k-1)+tn(k))/2))); else g(j)=v*f(xn(j+1),((tn(k-1)+tn(k))/2));

end end end for i=1:n-1 l(i)=u(k-1,i+1); end b=p*(l')+(g'); %解方程 ul=q\b; for i=2:n u(k,i)=ul(i-1); end end %得出误差 for i=1:m+1 for j=1:n+1 e1(i,j)=abs((f1(xn(j),tn(i)))-u(i,j)); end end e=e1'; for i=1:0.1/v:m+1 s((i-1)/(0.1/v)+1)=u(i,n/2+1); end s=s'; end %右端项 function f=f(x,t) f=-exp(x)*(cos(1/2-t)+2*sin(1/2-t)); end %精确解 function w=f1(x,t) w=exp(x)*sin(1/2-t); end