偏微分方程数值解大作业目录第一题 (3)第二题 (7)第三题 (16)第四题 (20)第五题 (26)第六题（附加题1） (39)第七题（附加题2） (45)第八题（附加题3） (51)第一题习题13.（1）解曲线图图1 （2）误差曲线图图2（3）表格表1 部分点处精确解和取不同步长时所得的数值解表2 取不同步长时部分结点处数值解的误差的绝对值和数值解的最大误差（4）MATLAB源代码M=64;a=0;b=pi/2;h=(b-a)/M;x=[a+h:h:b-h];u=zeros(M-1,M-1);u(1,1)=(2/h^2)+(x(1)-1/2)^2;u(1,2)=-(1/h^2);u(M-1,M-1)=(2/h^2)+(x(M-1)-1/2)^2;u(M-1,M-2)=-(1/h^2);for i=2:M-2u(i,i-1)=-(1/h^2);u(i,i)=(2/h^2)+(x(i)-1/2)^2;u(i,i+1)=-(1/h^2);endf=zeros(M-1,1)f(1)=(x(1).*x(1)-x(1)+5/4).*sin(x(1));f(M-1)=(x(M-1).*x(M-1)-x(M-1)+5/4).*sin(x(M-1))+1/h^2; for j=2:M-2f(j)=(x(j).*x(j)-x(j)+5/4).*sin(x(j));endy=inv(u)*f; true=sin(x); plot(x,y'-true)第二题习题二（P67 第3题）（1）h=1/4, τ=1/4精确解数值解误差（2）h=1/8, τ=1/8精确解数值解误差（3）h=1/16, τ=1/16精确解数值解误差（4）h=1/32, τ=1/32精确解数值解误差（5）h=1/64, τ=1/64精确解精确解误差（6）表格（7）Matlab代码function[p,e,u,x,y,k]=fivepoint(h,m,n,kmax,ep)%五点差分法和G-S迭代法解椭圆型方程%kmax为最大迭代次数；%m,n分别为x,y方向的网格数；%ep为精度；%u为差分解，p为精确解，e为误差；%例如在命令行窗口输入[p,e,u,x,y,k]=fivepoint(1/64,64,64,10000,1e-10);%代表步长为1/64,精度为10^(-10),最大迭代次数为10000的五点差分格式syms temp;u=zeros(n+1,m+1);x=0+(0:m)*h;y=0+(0:n)*h;for(j=1:n+1)u(j,1)=sin(y(j))+cos(y(j));u(j,m+1)=exp(1)*(sin(y(j))+cos(y(j)));endfor(i=1:m+1)u(1,i)=exp(x(i));u(n+1,i)=exp(x(i))*(sin(1)+cos(1));endt=zeros(n-1,m-1);for(k=1:kmax)for(j=2:n)for(i=2:m)temp=(u(j,i+1)+u(j,i-1)+u(j+1,i)+u(j-1,i))/4; t(j,i)=(temp-u(j,i))*(temp-u(j,i));u(j,i)=temp;endendt(j,i)=sqrt(t(j,i));if(k>kmax)break;endif(max(max(t))<ep)break;endendfor(j=1:n+1)for(i=1:m+1)p(j,i)=(sin(y(j))+cos(y(j)))*exp(x(i));e(j,i)=abs(u(j,i)-p(j,i));endend代码使用说明：在命令行窗口输入[p,e,u,x,y,k]=fivepoint(1/64,64,64,10000,1e-10);代表步长为1/64,精度为10^(-10),最大迭代次数为10000的五点差分格式surf(x,y,p)可作出近似解曲线surf(x,y,u)可作出精确解曲线surf(x,y,e)可作出误差曲线注意精度不要选取的太低，否则随着步长减小误差反而增大。

三张图分别为步长为1/4,1/16,1/64的误差曲线（精度均为10^(-10)）。

第三题习题三P138 第12题（1）h=0.010 tau=0.010k (x,t) 数值解精确解误差10 (0.5,0.1) 0.641417 0.642042 6.256223e-0420 (0.5,0.2) 0.486439 0.487230 7.914001e-0430 (0.5,0.3) 0.326667 0.327550 8.834544e-0440 (0.5,0.4) 0.163640 0.164597 9.575961e-0450 (0.5,0.5) -0.001021 0.000000 1.020916e-0360 (0.5,0.6) -0.165671 -0.164597 1.073862e-0370 (0.5,0.7) -0.328666 -0.327550 1.116054e-0380 (0.5,0.8) -0.488378 -0.487230 1.147092e-0390 (0.5,0.9) -0.643209 -0.642042 1.166668e-03100 (0.5,1.0) -0.791614 -0.790439 1.174587e-03h=0.005 tau=0.005k (x,t) 数值解精确解误差20 (0.5,0.1) 0.641729 0.642042 3.129321e-0440 (0.5,0.2) 0.486834 0.487230 3.959919e-0460 (0.5,0.3) 0.327108 0.327550 4.420366e-0480 (0.5,0.4) 0.164118 0.164597 4.790797e-04100 (0.5,0.5) -0.000511 0.000000 5.107000e-04120 (0.5,0.6) -0.165135 -0.164597 5.371294e-04140 (0.5,0.7) -0.328109 -0.327550 5.581797e-04160 (0.5,0.8) -0.487804 -0.487230 5.736512e-04180 (0.5,0.9) -0.642626 -0.642042 5.833907e-04200 (0.5,1.0) -0.791026 -0.790439 5.873012e-04h=0.0050 tau=0.0005k (x,t) 数值解精确解误差200 (0.5,0.1) 0.642011 0.642042 3.115887e-05400 (0.5,0.2) 0.487191 0.487230 3.948351e-05600 (0.5,0.3) 0.327506 0.327550 4.412446e-05800 (0.5,0.4) 0.164550 0.164597 4.786762e-051000 (0.5,0.5) -0.000051 0.000000 5.106902e-051200 (0.5,0.6) -0.164651 -0.164597 5.375135e-05 1400 (0.5,0.7) -0.327606 -0.327550 5.589538e-05 1600 (0.5,0.8) -0.487288 -0.487230 5.748076e-05 1800 (0.5,0.9) -0.642101 -0.642042 5.849178e-05 2000 (0.5,1.0) -0.790498 -0.790439 5.891837e-05h τ E∞(h, τ) E∞(2h,2τ)/E∞(h, τ)1/10 1/10 1.175210e-02 *1/20 1/20 5.910745e-03 1.9881/40 1/40 2.955702e-03 1.9991/80 1/80 1.478257e-03 1.9991/160 1/160 7.391688e-04 1.999 (2) t=1时数值解的误差曲线(3)外推步长1/100 步长1/200 数值解精确解误差（3）Matlab代码：h=input('输入空间步长h=');tau=input('输入时间步长tau=');%h=1/100;tau=1/200;r=tau/h^2;m=1/h; n=1/tau;i=1:m-1; k=1:n-1;A=diag(-r*ones(m-2,1),-1)+diag((1+2*r)*ones(m-1,1))+diag(-r*ones(m-2,1),1);B=diag(r*ones(m-2,1),-1)+diag((1-2*r)*ones(m-1,1))+diag(r*ones(m-2,1),1);%C=diag(-1/2*s^2*ones(m-2,1),-1)+diag((2+s^2)*ones(m-1,1))+diag(-1/2*s^2*ones(m-2,1),1);xi=0+1/m:1/m:1-1/m;tk=1/n:1/n:1;ui0=(exp(xi)*sin(1/2))'; u00=sin(1/2); um0=exp(1)*sin(1/2);u0k=sin(1/2-tk); umk=exp(1)*sin(1/2-tk);fik=zeros(m-1,n);for j=1:nfik(:,j)=-exp(xi)*(cos(1/2-tk(j))+2*sin(1/2-tk(j)));enduik=zeros(m-1,n);uik(:,1)=inv(A)*(B*ui0+[r*u00;zeros(m-3,1);r*um0]+[r*u0k(1);zeros(m-3,1);r*umk(1)]+tau*(-exp(xi)*(cos(1/2)+2*sin(1/2)))');%first k=1%uik(:,2)=inv(A)*(B*uik(:,1)+[r*u0k(1);zeros(m-3,1);r*umk(1)]+[r*u0k(2);zeros(m-3,1);r*umk(2)]+tau*fik(:,1));%%%%%%%%%%%%%%%%for k=2:n-1D=(B*uik(:,k)+[r*u0k(k);zeros(m-3,1);r*umk(k)]+[r*u0k(k+1);zeros(m-3,1);r*umk(k+1)]+tau*fik(:,k));uik(:,k+1)=inv(A)*D;endtrue=zeros(m-1,n);for i=1:m-1true(i,:)=exp(xi(i))*sin(1/2-tk);endfprintf('h=%.4f tau=%.4f\n',h,tau);fprintf('k (x,t) 数值解精确解误差\n');for k=1:10p=int16(k*0.1/tau);if p<100fprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));elsefprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));endenda=max(abs(uik-true));E=max(a);fprintf('Emax=%e',E);[x,t]=meshgrid(xi,tk);mesh(x,t,abs(uik-true)');xlabel('xi');ylabel('tk');hold on;%e=max(true-uik);% h=max(e);% plot(e)%surf(abs(uik-true))%surf(true)% plot(xi,abs(uik(:,n)-true(:,n)))% hold on;第四题P176 习题四第8题(1)表格h=0.010 tau=0.005(x,t) 数值解精确解误差(0.5,0.1) 0.049979 0.049979 4.791843e-08 (0.5,0.2) 0.099833 0.099833 7.083966e-08 (0.5,0.3) 0.149438 0.149438 4.406409e-08 (0.5,0.4) 0.198669 0.198669 5.641112e-08 (0.5,0.5) 0.247404 0.247404 2.789924e-07 (0.5,0.6) 0.295519 0.295520 9.784169e-07 (0.5,0.7) 0.342896 0.342898 1.775191e-06 (0.5,0.8) 0.389416 0.389418 2.623113e-06 (0.5,0.9) 0.434962 0.434966 3.474007e-06 (0.5,1.0) 0.479421 0.479426 4.279944e-06h=0.005 tau=0.010(x,t) 数值解精确解误差(0.5,0.1) 0.049979 0.049979 1.917987e-07(0.5,0.2) 0.099834 0.099833 2.836035e-07 (0.5,0.3) 0.149438 0.149438 1.765491e-07 (0.5,0.4) 0.198669 0.198669 2.261917e-07 (0.5,0.5) 0.247403 0.247404 1.148338e-06 (0.5,0.6) 0.295516 0.295520 3.962614e-06 (0.5,0.7) 0.342891 0.342898 7.167473e-06 (0.5,0.8) 0.389408 0.389418 1.059995e-05 (0.5,0.9) 0.434951 0.434966 1.405730e-05 (0.5,1.0) 0.479408 0.479426 1.748442e-05h τE∞(h, τ) E∞(2h,2τ)/E∞(h, τ)1/10 1/10 1.767294e-03 *1/20 1/20 4.356573e-04 4.057 1/40 1/40 1.100309e-04 3.959 1/80 1/80 2.753878e-05 3.995 1/160 1/160 6.882379e-06 4.001 1/320 1/320 1.720007e-06 4.001(2)精确解 (h=1/100, τ=1/100)数值解(h=1/10, τ=1/10)(3)上曲面 (h=1/10, τ=1/10) 中曲面(h=1/20, τ=1/20) 下曲面(h=1/40, τ=1/40)（4）Matlab代码：h=input('输入空间步长h=');tau=input('输入时间步长tau='); s=tau/h;m=1/h; n=1/tau;i=1:m-1; k=1:n-1;A=diag(-1/2*s^2*ones(m-2,1),-1)+diag((1+s^2)*ones(m-1,1))+diag(-1/2*s^2*ones(m-2,1),1);B=diag(1/2*s^2*ones(m-2,1),-1)+diag(-(1+s^2)*ones(m-1,1))+diag(1/2*s^2*ones(m-2,1),1);C=diag(-1/2*s^2*ones(m-2,1),-1)+diag((2+s^2)*ones(m-1,1))+diag(-1/2*s^2*ones(m-2,1),1);xi=0+1/m:1/m:1-1/m;tk=1/n:1/n:1;ui0=zeros(m-1,1); u00=0; um0=0;u0k=zeros(n,1); umk=sin(tk);fik=zeros(m-1,n);for j=1:nfik(:,j)=(((tk(j)).^2-(xi).^2).*sin((xi).*(tk(j))))'; enduik=zeros(m-1,n);uik(:,1)=inv(C)*(2*tau*xi'+[zeros(m-2,1);1/2*s^2*umk(1)]);%first k=1%D=2*uik(:,1)+B*ui0+tau^2*fik(:,1)+[zeros(m-2,1);1/2*s^2*umk(2)];uik(:,2)=inv(A)*D;%%%%%%%%%%%%%%%%for k=2:n-1D=2*uik(:,k)+B*uik(:,k-1)+[zeros(m-2,1);1/2*s^2*umk(k-1)]+tau^2*fik(:,k)+[zeros(m-2,1);1/2*s^2*umk(k+1)];uik(:,k+1)=inv(A)*D;endtrue=zeros(m-1,n);for i=1:m-1true(i,:)=sin(xi(i)*tk);endfprintf('h=%.3f tau=%.3f\n',h,tau);fprintf('k (x,t) 数值解精确解误差\n');for k=1:10p=int16(k*0.1/tau);if p<100fprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),tru e(m/2,p),abs(uik(m/2,p)-true(m/2,p)));elsefprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),tru e(m/2,p),abs(uik(m/2,p)-true(m/2,p)));endenda=max(abs(uik-true));E=max(a);fprintf('Emax=%e',E);[x,t]=meshgrid(xi,tk);mesh(x,t,true');xlabel('xi');ylabel('tk');hold on;%e=max(true-uik);% h=max(e)% plot(e)%surf(abs(uik-true))%surf(true)%plot(abs(uik(:,n)-true(:,n)))第五题习题56.（1）曲线图① h=0.05 τ= 0.05数值解精确解误差曲线图② h=0.025 τ= 0.025数值解精确解误差曲线图（2）表格① h=0.0500 τ=0.0500k (x,y,t) 数值解精确解误差02 (0.5,0.5,0.1) 0.022374 0.024997 2.623812e-03 04 (0.5,0.5,0.2) 0.049668 0.049979 3.116082e-04 06 (0.5,0.5,0.3) 0.074901 0.074930 2.839556e-05 08 (0.5,0.5,0.4) 0.099852 0.099833 1.849222e-05 10 (0.5,0.5,0.5) 0.124713 0.124675 3.871957e-05 12 (0.5,0.5,0.6) 0.149497 0.149438 5.852752e-05 14 (0.5,0.5,0.7) 0.174189 0.174108 8.110513e-05 16 (0.5,0.5,0.8) 0.198776 0.198669 1.066039e-04 18 (0.5,0.5,0.9) 0.223241 0.223106 1.347368e-04 20 (0.5,0.5,1.0) 0.247569 0.247404 1.651274e-04② h=0.0250 τ=0.0250k (x,y,t) 数值解精确解误差04 (0.5,0.5,0.1) 0.024101 0.024997 8.967699e-04 08 (0.5,0.5,0.2) 0.049855 0.049979 1.240867e-04 12 (0.5,0.5,0.3) 0.074916 0.074930 1.412105e-05 16 (0.5,0.5,0.4) 0.099837 0.099833 3.688824e-06 20 (0.5,0.5,0.5) 0.124684 0.124675 9.652543e-06 24 (0.5,0.5,0.6) 0.149453 0.149438 1.475416e-05 28 (0.5,0.5,0.7) 0.174129 0.174108 2.044073e-05 32 (0.5,0.5,0.8) 0.198696 0.198669 2.684217e-05 36 (0.5,0.5,0.9) 0.223140 0.223106 3.389905e-05 40 (0.5,0.5,1.0) 0.247445 0.247404 4.151868e-05（3）MATLAB源代码clear;%%h=input('输入空间步长h=');tau=input('输入时间步长tau=');r=tau/h^2;m=1/h;n=1/tau;i=1:m-1;j=1:m-1;k=1:n-1;A=diag(-r/2*ones(m-2,1),-1)+diag((1+r)*ones(m-1,1))+diag(-r/2*ones(m-2,1),1);B=diag(r/2*ones(m-2,1),-1)+diag((1-r)*ones(m-1,1))+diag(r/2*ones(m-2,1),1); xi=0+1/m:1/m:1-1/m;yj=0+1/m:1/m:1-1/m;tk=1/n:1/n:1;uij0=zeros(m-1,m-1);u000=0;u0m0=0;um00=0;umm0=0;uim0=zeros(m-1,1);um0k=0;ummk=sin(tk);umj0=0;ui00=zeros(m-1,1);%%u0jk=zeros(m-1,n);umjk=zeros(m-1,n);ui0k=zeros(m-1,n);uimk=zeros(m-1,n);for k=1:nu0jk(:,k)=0;umjk(:,k)=sin(yj*tk(k))';ui0k(:,k)=0;uimk(:,k)=sin(xi*tk(k))';endfijk=zeros(m-1,m-1,n);fij0=0;for k=1:nfor j=1:m-1fijk(:,j,k)=((xi.^2+yj(j).^2).*(tk(k).^2).*sin(xi*yj(j)*tk(k))+xi.*yj(j).*cos( xi*yj(j)*tk(k)))';endenduijk0p5=zeros(m-1,m-1,n);uijk=zeros(m-1,m-1,n);true=zeros(m-1,m-1,n);for k=1:nfor j=1:m-1true(:,j,k)=sin(xi*yj(j)*tk(k))';endend%%for k=1:nfor j=1:m-1if j==1if k==1u0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*0))-1/4*tau*1/h^2*((um0k-2*umjk(j,k)+umjk(j+1,k)-(um0k-2*umj0+umj0)));elseu0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((um0k-2*umjk(j,k)+umjk(j+1,k)-(um0k-2*umjk(j,k-1)+umjk(j+1,k-1))));endelseif j==m-1if k==1u0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+0)-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+ummk(k)-(umj0-2*umj0+umm0)));elseu0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+ummk(k)-(umjk(j-1,k-1)-2*umjk(j,k-1)+ummk(k-1))));endelseif k==1u0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+0)-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+umjk(j+1,k)-(umj0-2*umj0+umj0)));elseu0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+umjk(j+1,k)-(umjk(j-1,k-1)-2*umjk(j,k-1)+umjk(j+1,k-1))));endendif k==1if j==1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uij0(p,1);if p~=1temp1(p,p-1)=uij0(p-1,2);endendtemp2=ui00+[zeros(m-2,1);uij0(m-1,2)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%%%%%%%%%%%%%%%%%%%%%addtemp3=(fijk(:,j,k)+fij0)/2;temp4=diag((1-r)*ones(m-1,1))+diag(r/2*ones(m-2,1),1);uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3); elseif j==m-1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uij0(p,m-2);if p~=1temp1(p,p-1)=uij0(p-1,m-1);endendtemp2=uim0+[zeros(m-2,1);2*(1-r)/r*uij0(m-1,m-1)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%%%%%%%%%%%%%%%%%%%%%%%addtemp3=(fijk(:,j,k)+fij0)/2;temp4=diag(r/2*ones(m-1,1))+diag((1-r)*ones(m-2,1),1);uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3); elsetemp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uij0(p,j);if p~=1temp1(p,p-1)=uij0(p-1,j+1);endif p~=m-1temp1(p,p+1)=uij0(p+1,j-1);endendtemp2=[uij0(1,j-1);zeros(m-3,1);uij0(m-1,j+1)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%%%%%%%%%%%%%%%%%%%addtemp3=(fijk(:,j,k)+fij0)/2;temp4=B;uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3); endelseif j==1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk(p,1,k-1);if p~=1temp1(p,p-1)=uijk(p-1,2,k-1);endendtemp2=ui0k(:,k-1)+[zeros(m-2,1);uijk(m-1,2,k-1)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%%%%%%%%%%%%%%addtemp3=(fijk(:,j,k)+fijk(:,j,k-1))/2;temp4=diag((1-r)*ones(m-1,1))+diag(r/2*ones(m-2,1),1);uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3);elseif j==m-1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk(p,m-2,k-1);if p~=1temp1(p,p-1)=uijk(p-1,m-1,k-1);endendtemp2=uimk(:,k-1)+[zeros(m-2,1);2*(1-r)/r*uijk(m-1,m-1,k-1)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%all add in temp2temp3=(fijk(:,j,k)+fijk(:,j,k-1))/2;temp4=diag(r/2*ones(m-1,1))+diag((1-r)*ones(m-2,1),1);uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3); elsetemp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk(p,j,k-1);if p~=1temp1(p,p-1)=uijk(p-1,j+1,k-1);endif p~=m-1temp1(p,p+1)=uijk(p+1,j-1,k-1);endendtemp2=[uijk(1,j-1,k-1);zeros(m-3,1);uijk(m-1,j+1,k-1)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%%%%%%%%%%%%%addtemp3=(fijk(:,j,k)+fijk(:,j,k-1))/2;temp4=B;uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3); endendend%%%%%%%%%%%%%%%%%%%%%%edge Vec%%%%%%%%%%%%%%%%%%%u0jk0p5vec=0;umjk0p5vec=zeros(m-1,n);%%used before and now is freeif k==1for j=1:m-1if j==1u0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+0)-1/4*tau*1/h^2*((um0k-2*umjk(j,k)+umjk(j+1,k)-(umj0-2*umj0+umj0)));elseif j==m-1u0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+0)-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+ummk(k)-(umj0-2*umj0+umj0)));elseu0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+0)-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+umjk(j+1,k)-(umj0-2*umj0+umj0))); endendelsefor j=1:m-1if j==1u0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((um0k-2*umjk(j,k)+umjk(j+1,k)-(um0k-2*umjk(j,k-1)+umjk(j+1,k-1))));elseif j==m-1u0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+ummk(k)-(umjk(j-1,k-1)-2*umjk(j,k-1)+ummk(k-1))));elseu0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+umjk(j+1,k)-(umjk(j-1,k-1)-2*umjk(j,k-1)+umjk(j+1,k-1))));endendend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for i=1:m-1if i==1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk0p5(1,p,k);if p~=1temp1(p,p-1)=uijk0p5(2,p-1,k);endendtemp2=u0jk0p5vec+[zeros(m-2,1);uijk0p5(2,m-1,k)]+[ui0k(1,k);zeros(m-3,1);uimk(1,k)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if k==1temp3=(fijk(i,:,k)+fij0)'/2;elsetemp3=(fijk(i,:,k)+fijk(i,:,k-1))'/2;endtemp4=diag((1-r)*ones(m-1,1))+diag(r/2*ones(m-2,1),1);uijk(i,:,k)=(A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3))';elseif i==m-1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk0p5(m-2,p,k);if p~=1temp1(p,p-1)=uijk0p5(m-1,p-1,k);endendtemp2=umjk0p5vec(:,k)+[zeros(m-2,1);(1-r)*2/r*uijk0p5(m-1,m-1,k)]+[ui0k(m-1,k);zeros(m-3,1);uimk(m-1,k)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if k==1temp3=(fijk(i,:,k)+fij0)'/2;elsetemp3=(fijk(i,:,k)+fijk(i,:,k-1))'/2;endtemp4=diag((r/2)*ones(m-1,1))+diag((1-r)*ones(m-2,1),1);uijk(i,:,k)=(A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3))';elsetemp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk0p5(i,p,k);if p~=1temp1(p,p-1)=uijk0p5(i+1,p-1,k);endif p~=m-1temp1(p,p+1)=uijk0p5(i-1,p+1,k);endendtemp2=[uijk0p5(i-1,1,k);zeros(m-3,1);uijk0p5(i+1,m-1,k)]+[ui0k(i,k);zeros(m-3,1);uimk(i,k)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if k==1temp3=(fijk(i,:,k)+fij0)'/2;elsetemp3=(fijk(i,:,k)+fijk(i,:,k-1))'/2;endtemp4=B;uijk(i,:,k)=(A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3))';endendendsurf(uijk(:,:,n));figure(2)surf(true(:,:,n));figure(3)surf(uijk(:,:,n)-true(:,:,n));fprintf('h=%.4f tau=%.4f\n',h,tau);fprintf('k (x,y,t) 数值解精确解误差 \n');for k=1:10p=int16(k*0.1/tau);if p<100fprintf('%d(0.5,0.5,%.1f) %f %f %e\n',p,0.1*k,uijk(m/2,m/2,p),true(m/2,m/2 ,p),abs(uijk(m/2,m/2,p)-true(m/2,m/2,p)));elsefprintf('%d(0.5,0.5,%.1f) %f %f %e\n',p,0.1*k,uijk(m/2,m/2,p),true(m/2,m/2 ,p),abs(uijk(m/2,m/2,p)-true(m/2,m/2,p)));endend第六题补充题1（教材P94/算例3.4）（1）差分格式ðu ðt (x i,t k)−að2uðt2(x i,t k)=f(x i,t k)u i k+1+u i k−12τ−aδx2u i k−1+δx2u i k+12=f(x i,t k)（2）算例3.4部分结点数值解精确解和误差的绝对值① h=0.100 τ=0.010k (x,t) 数值解精确解误差10 (0.5,0.1) 1.822237 1.822119 1.181012e-0420 (0.5,0.2) 2.013930 2.013753 1.773622e-0430 (0.5,0.3) 2.225754 2.225541 2.135582e-0440 (0.5,0.4) 2.459846 2.459603 2.425890e-0450 (0.5,0.5) 2.718552 2.718282 2.705635e-0460 (0.5,0.6) 3.004466 3.004166 2.999407e-0470 (0.5,0.7) 3.320449 3.320117 3.318310e-0480 (0.5,0.8) 3.669664 3.669297 3.668593e-0490 (0.5,0.9) 4.055605 4.055200 4.054907e-04100 (0.5,1.0) 4.482137 4.481689 4.481546e-04② h=0.050 τ=1/400k (x,t) 数值解精确解误差40 (0.5,0.1) 1.822148 1.822119 2.874976e-05 80 (0.5,0.2) 2.013796 2.013753 4.313840e-05 120 (0.5,0.3) 2.225593 2.225541 5.191899e-05 160 (0.5,0.4) 2.459662 2.459603 5.896391e-05 200 (0.5,0.5) 2.718348 2.718282 6.575685e-05 240 (0.5,0.6) 3.004239 3.004166 7.289348e-05 280 (0.5,0.7) 3.320198 3.320117 8.064224e-05 320 (0.5,0.8) 3.669386 3.669297 8.915426e-05 360 (0.5,0.9) 4.055299 4.055200 9.854220e-05 400 (0.5,1.0) 4.481798 4.481689 1.089103e-04③取不同步长时数值解的最大误差④误差曲线图⑤误差曲面图图中：上曲面h=1/10 τ=1/100下曲面h=1/20 τ=1/400（3）MATLAB源代码h=input('输入空间步长h=');tau=input('输入时间步长tau=');r=tau/h^2;m=1/h;n=1/tau;i=1:m-1;k=1:n-1;A=diag(-r*ones(m-2,1),-1)+diag((1+2*r)*ones(m-1,1))+diag(-r*ones(m-2,1),1);B=diag(r*ones(m-2,1),-1)+diag((1-2*r)*ones(m-1,1))+diag(r*ones(m-2,1),1);xi=0+1/m:1/m:1-1/m;tk=1/n:1/n:1;ui0=exp(xi)';u00=1;%u01=1+tau;um0=exp(1);u0k=exp(tk)';umk=exp(1+tk)';fik=zeros(m-1,1);uik=zeros(m-1,n);%uik2=ones(m-1,n-1);uik(:,1)=(1+tau)*exp(xi);%first k=1%C=B*ui0+2*tau*fik+[r*(u0k(2)+r*u00);zeros(m-3,1);r*(umk(2)+um0)];%uik(:,2)=tridiag(A,C,m-1);uik(:,2)=inv(A)*C;%%%%%%%%%%%%%%%%for k=2:n-1C=B*uik(:,k-1)+2*tau*fik+[r*(u0k(k+1)+u0k(k-1));zeros(m-3,1);r*(umk(k+1)+umk(k-1))];uik(:,k+1)=inv(A)*C;endU2=exp(xi+1)';true=zeros(m-1,n);for i=1:m-1true(i,:)=exp(xi(i)+tk);endfprintf('h=%.3f tau=%.3f\n',h,tau);fprintf('k (x,t) 数值解精确解误差\n');for k=1:10p=int16(k*0.1/tau);if p<100fprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));elsefprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));endenda=max(abs(uik-true));E=max(a);fprintf('Emax=%e',E);% e=max(true-uik); % h=max(e)%plot(e)%plot(abs(uik(:,n)-true(:,n)))第七题补充题2 非线性方程（教材P129/算例3.4）（1）差分格式ðu ðt −u4ð2uðt2=f(x,t)u i k+1+u i k−12τ−(u i k)4δx2(u i k−1+u i k+12)=f(x i,t k)（2）取不同步长时数值解的最大误差（3）误差曲线图（4）误差曲面图图中：上曲面h=τ=1/10中曲面h=τ=1/20下曲面h=τ=1/40（5）MATLAB源代码h=input('输入空间步长h=');tau=input('输入时间步长tau=');r=tau/h^2;m=1/h;n=1/tau;i=1:m-1;k=1:n-1;xi=0+1/m:1/m:1-1/m;tk=1/n:1/n:1;ui0=sin(pi/4+xi*pi/2)';u00=sqrt(2)/2;um0=sqrt(2)/2;u0k=sqrt(2)/2*exp(-tk)';umk=sqrt(2)/2*exp(-tk)';fik=zeros(m-1,n);for j=1:nfik(:,j)=(exp(-tk(j))*sin(pi/2+xi*pi/2).*(pi^2/4*exp(-4*tk(j)).*sin(pi/4+xi*pi/2).^4-1))'; enduik=zeros(m-1,n);%uik2=ones(m-1,n-1);uik(:,1)=(sin(pi/4+xi*pi/2)+tau*(sin(pi/4+xi*pi/2).^4).*(-pi^2/4*sin(pi/4+xi*pi/2))+tau*sin(pi/4+xi*pi/2).*(pi^2/4*sin(pi/4+xi*pi/2).^4-1))';%first k=1%A=zeros(m-1,m-1);B=zeros(m-1,m-1);A(1,1)=1+2*r*(uik(1,1)^4);A(1,2)=-r*(uik(1,1)^4);B(1,1)=1-2*r*(uik(1,1)^4);B(1,2)=r*(uik(1,1)^4);A(m-1,m-2)=-r*(uik(m-1,1)^4);A(m-1,m-1)=1+2*r*(uik(m-1,1)^4);B(m-1,m-2)=r*(uik(m-1,1)^4);B(m-1,m-1)=1-2*r*(uik(m-1,1)^4);for i=2:m-2A(i,i-1)=-r*(uik(i,1)^4);A(i,i)=1+2*r*(uik(i,1)^4);A(i,i+1)=-r*(uik(i,1)^4);B(i,i-1)=r*(uik(i,1)^4);B(i,i)=1-2*r*(uik(i,1)^4);B(i,i+1)=r*(uik(i,1)^4);endC=B*ui0+2*tau*fik(:,1)+[r*(uik(1,1)^4*(u0k(2)+u00));zeros(m-3,1);r*(uik(m-1,1))^4*(umk(2)+um0)];%uik(:,2)=tridiag(A,C,m-1);uik(:,2)=inv(A)*C;%%%%%%%%%%%%%%%%for k=2:n-1A=zeros(m-1,m-1);B=zeros(m-1,m-1);A(1,1)=1+2*r*(uik(1,k)^4);A(1,2)=-r*(uik(1,k)^4);B(1,1)=1-2*r*(uik(1,k)^4);B(1,2)=r*(uik(1,k)^4);A(m-1,m-2)=-r*(uik(m-1,k)^4);A(m-1,m-1)=1+2*r*(uik(m-1,k)^4);B(m-1,m-2)=r*(uik(m-1,k)^4);B(m-1,m-1)=1-2*r*(uik(m-1,k)^4);for i=2:m-2A(i,i-1)=-r*(uik(i,k)^4);A(i,i)=1+2*r*(uik(i,k)^4);A(i,i+1)=-r*(uik(i,k)^4);B(i,i-1)=r*(uik(i,k)^4);B(i,i)=1-2*r*(uik(i,k)^4);B(i,i+1)=r*(uik(i,k)^4);endC=B*uik(:,k-1)+2*tau*fik(:,k)+[r*(uik(1,k)^4)*(u0k(k+1)+u0k(k-1));zeros(m-3,1);r*(uik(m-1,k)^4)*(umk(k+1)+umk(k-1))];uik(:,k+1)=inv(A)*C;endtrue=zeros(m-1,n);for i=1:m-1true(i,:)=exp(-tk).*(sin(pi/4+xi(i)*pi/2));end% e=max(true-uik);% h=max(e)% plot(e)% plot(uik(:,90));% hold on;plot(true(:,end)-uik(:,end))%surf(abs(uik-true))fprintf('h=%.3f tau=%.3f\n',h,tau);fprintf('k (x,t) 数值解精确解误差\n'); for k=1:10p=int16(k*0.1/tau);if p<100fprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));elsefprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));endenda=max(abs(uik-true));E=max(a);fprintf('Emax=%e',E);。