偏微分方程程序附件.doc

程序一、 计算下面双曲方程的近似解解法：1、 一阶迎风格式程序：function upwind1(h,lamda) %迎风格式 一阶精度 显示格式x=h:h:10-h;tao=lamda*h;t=tao:tao:5-tao;X,T=meshgrid(x,t);i=1:10/h;j=1:5/tao;u(i,j)=0;for n=2:5/tao for m=2:10/h u(m,n)=u(m,n-1)-lamda*(u(m,n-1)-u(m-1,n-1)+2*tao*(n-1)*tao*sin(m*h)+(n-1)2*tao3*cos(m*h); %迎风格式 endendfor n=1:5/tao-1 for m=1:10/h-1 U(m,n)=u(m+1,n+1); endendU=U;figure(1);grid on;mesh(X,T,U);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(迎风格式差分曲面图);%计算误差ii=1:10/h-1;jj=1:5/tao-1;uu(ii,jj)=0;for nn=1:5/tao-1 for mm=1:10/h-1 uu(mm,nn)=(nn*tao)2*sin(mm*h); endendUU=uu;error= U-UU;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(迎风格式误差曲面图);figure(3);mesh(X,T,UU);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(古典解);2、 Lax-Wendroff格式function lax_w1(h,lamda) %Lax-Wendroff 二阶精度 显示格式x=h:h:10-h;tao=lamda*h;t=tao:tao:5-tao;X,T=meshgrid(x,t);i=1:10/h;j=1:5/tao;u(i,j)=0;for n=2:5/tao for m=2:10/h if m=10/h u(m,n)=u(m,n-1)-lamda*(u(m,n-1)-u(m-1,n-1)+2*tao*(n-1)*tao*sin(m*h). +(n-1)2*tao3*cos(m*h); %用迎风格式处理边界问题 else u(m,n)=u(m,n-1)-0.5*lamda*(u(m+1,n-1)-u(m-1,n-1)+2*tao*(n-1)*tao*sin(m*h). +tao*(n-1)*tao)2*cos(m*h)+0.5*lamda2*(u(m+1,n-1)-2*u(m,n-1)+u(m-1,n-1).+tao2*(n-1)*tao*(sin(m*h)-cos(m*h)+0.5*tao2*(n-1)2*tao2*(sin(m*h)+cos(m*h); %Lax-Wendroff差分格式 end endendfor n=1:5/tao-1 for m=1:10/h-1 U(m,n)=u(m+1,n+1); endendU=U;figure(1);grid on;mesh(X,T,U);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(Lax-Wendroff格式差分曲面图);%计算误差ii=1:10/h-1;jj=1:5/tao-1;uu(ii,jj)=0;for nn=1:5/tao-1 for mm=1:10/h-1 uu(mm,nn)=(nn*tao)2*sin(mm*h); endendUU=uu;error=U-UU;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(Lax-Wendroff格式误差曲面图);3、 隐式中心格式function implicit_1(h,lamda) %隐式格式，精度为：时间一阶，空间二阶；无条件稳定x=h:h:10-h;tao=lamda*h;t=tao:tao:5-tao;X,T=meshgrid(x,t);time1=5/tao; %时间为1时的时间层网格点space1=10/h; %空间为1时的时间层网格点i=1:space1;j=1:time1;u(i,j)=0;r=lamda/2;%使用追赶法求每一时间层的节点值k=1:space1-1; %初始化系数矩阵 a(k,k)=0; a(1,1)=1; a(1,2)=r;for k=2:space1-2 %定义系数矩阵 a(k,k-1)=-r; a(k,k)=1; a(k,k+1)=r;end a(space1-1,space1-1)=1; p=1:space1-1; q=1:time1-1; Z(p,q)=0; %矩阵a*x(i)=Z(i); for n=2:time1 for m=2:space1-1 Z(m-1,n-1)=u(m,n-1)+tao*(2*(n-1)*tao*sin(m-1)*h)+. n2*tao2*cos(m*h); end Z(space1-1,n-1)=u(space1,n-1)-lamda*(u(space1,n-1)-u(space1-1,n-1). +2*tao*(n-1)*tao*sin(space1*h)+. (n-1)2*tao3*cos(space1*h); %边界离散，时间一阶，空间一阶 E=chase(a,Z(:,n-1); %调用追赶法程序，程序在下面 for v=1:space1-1 u(v+1,n)=E(v); endend for n=1:5/tao-1 for m=1:10/h-1 U(m,n)=u(m+1,n+1); endendU=U;figure(1);grid on;mesh(X,T,U);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(隐式格式差分曲面图);%计算误差ii=1:10/h-1;jj=1:5/tao-1;uu(ii,jj)=0;for nn=1:5/tao-1 for mm=1:10/h-1 uu(mm,nn)=(nn*tao)2*sin(mm*h); endendUU=uu;error=UU-U;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(隐式格式误差曲面图);调用追赶法接三对角矩阵函数chase.m为：function X=chase(A,d)%追赶法只用来解三对角方程组，其中A为方程的系数，d为右端项%a为-1对角线上元素%b为主对角线上的元素%c为+1对角线上的元素a=diag(A,-1);b=diag(A);c=diag(A,1);n=length(b);U(1)=b(1);y(1)=d(1);for i=2:n L(i-1)=a(i-1)/U(i-1); U(i)=b(i)-c(i-1)*L(i-1); y(i)=d(i)-L(i-1)*y(i-1);endX(n)=y(n)/U(n);for i=n-1:-1:1 X(i)=(y(i)-c(i)*X(i+1)/U(i);end二、 计算下面双曲方程的近似解解法：1、 二阶显示格式程序：function explicit_2(h,lamda) %波动方程显示格式，二阶精度，lamda1时稳定tao=h*lamda;t=tao:tao:1-tao;x=h:h:1-h;time1=1/tao+1; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点X,T=meshgrid(x,t);i=1:space1;j=1:time1;u(i,j)=0;for n=2:time1 %边界条件离散 u(space1,n)=sin(n-1)*tao);endfor m=2:1/h %初始条件离散 二阶精度 u(m,2)=u(m,1)+tao*h*(m-1)+0.5*lamda2*(u(m+1,1)-2*u(m,1)+u(m-1,1);end for n=3:1/tao for m=2:space1-1 u(m,n)=2*u(m,n-1)-u(m,n-2)+lamda2*(u(m+1,n-1)-2*u(m,n-1)+u(m-1,n-1)+. tao2*(n-1)*tao)2-(m-1)*h)2)*sin(n-1)*tao*(m-1)*h); %波动方程显示格式 endendfor n=1:time1-2 for m=1:space1-2 U(m,n)=u(m+1,n+1); endendU=U;figure(1);grid on;mesh(X,T,U);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(波动方程显示格式差分曲面图);%计算古典解ii=1:space1-2;jj=1:time1-2;uu(ii,jj)=0;for nn=1:time1-2 for mm=1:space1-2 uu(mm,nn)=sin(mm*h*nn*tao); endendUU=uu;error= U-UU;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(波动方程显示格式误差曲面图);figure(3);mesh(X,T,UU);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(波动方程古典解);2、 隐式二阶格式程序：function implicit_2(h,lamda) %波动方程隐式格式，二阶精度，无条件稳定tao=h*lamda;t=tao:tao:1-tao;x=h:h:1-h;time1=1/tao+1; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点X,T=meshgrid(x,t);i=1:space1;j=1:time1;u(i,j)=0; %离散化的网格点v(i,j)=0; %定义关于时间的一阶导数for i=1:space1; v(i,1)=(i-1)*h; %时间的一阶导数初始值离散化endfor i=1:space1-1 for j=1:time1 w(i,j)=0; %定义关于空间的一阶导数 endend for i=1:space1-1 w(i,1)=0; %空间的一阶导数初始值离散化endfor n=2:time1 %边界条件离散 u(space1,n)=sin(n-1)*tao);end%使用追赶法求每一时间层v的节点值,需要定义系数矩阵a和非其次项Zr=-0.25*lamda2;m=1+0.5*lamda2;k=1:space1-2; %初始化系数矩阵 a(k,k)=0; a(1,1)=m; a(1,2)=r;for k=2:space1-3 %定义系数矩阵 a(k,k-1)=r; a(k,k)=m; a(k,k+1)=r;end a(space1-2,space1-3)=r; a(space1-2,space1-2)=m;%定义非齐次项 p=1:space1-2; q=1:time1-1; Z(p,q)=0; %矩阵a*x(i)=Z(i); %使用隐式格式计算 for jj=2:time1 %从第二时间层开始计算 v(1,jj)=0; %v左边界值离散 v(space1,jj)=(u(space1,jj)-u(space1,jj-1)/tao; %v右边界值离散 %计算非齐次项 for m=1:space1-3 Z(m,jj-1)=v(m+1,jj-1)+lamda*(w(m+1,jj-1)-w(m,jj-1)-. r*(v(m+2,jj-1)-2*v(m+1,jj-1)+v(m,jj-1); end Z(space1-2,jj-1)=-r*v(space1,jj)+v(space1-1,jj-1)+lamda*. (w(space1-1,jj-1)-w(space1-2,jj-1)-. r*(v(space1,jj-1)-2*v(space1-1,jj-1)+v(space1-2,jj-1); E=chase(a,Z(:,jj-1); for kk=1:space1-2 v(kk+1,jj)=E(kk); end %计算w(:,jj)和u(:,jj) for ii=1:space1-1 w(ii,jj)=0.5*lamda*(v(ii+1,jj)-v(ii,jj)+w(ii,jj-1)+. 0.5*lamda*(v(ii+1,jj-1)-v(ii,jj-1); end for g=2:space1-1 u(g,jj)=u(g,jj-1)+tao*v(g,jj); endendfor n=1:time1-2 for m=1:space1-2 U(m,n)=u(m+1,n+1); endendU=U;figure(1);grid on;mesh(X,T,U);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(波动方程显示格式差分曲面图);%计算古典解i1=1:space1-2;jj=1:time1-2;uu(i1,jj)=0;for nn=1:time1-2 for mm=1:space1-2 uu(mm,nn)=sin(mm*h*nn*tao); endendUU=uu;error= U-UU;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(波动方程显示格式误差曲面图);三、 计算下面抛物方程的近似解1、 向前差分格式程序：function diffusion_1(h,lamda,time) %扩散方程向前差分格式，时间一阶精度，空间二阶精度，lamda=0,5稳定%h为空间步长，lamda为网格比，t为演化时间%h=0.1;time=1;lamda=0.5;tao=h2*lamda;t=0:tao:time;x=0:h:1;X,T=meshgrid(x,t);time1=time/tao+2; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点i=1:space1;j=1:time1;u(i,j)=0; %离散化的网格点wucha(1,i)=0;for i=1:space1-1 %边界条件初始化 u(i,1)=sin(pi*(i-1)*h);endfor n=2:time1 for m=2:space1-1 u(m,n)=(1-2*lamda)*u(m,n-1)+lamda*u(m-1,n-1)+lamda*u(m+1,n-1); endendu=u;figure(1);grid on;mesh(X,T,u);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(向前差分格式曲面图);%计算误差ii=1:space1;jj=1:time1;uu=exp(-pi2.*T).*sin(pi.*X);error=uu-u;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(向前差分格式误差曲面图);%plot(error);figure(3);mesh(X,T,uu);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(古典解);wucha=error(22,:)figure(4)plot(wucha,ro-);xlabel(t=0.1s，x轴);ylabel(误差轴);title(0.1秒时差分格式解与精确解的误差);2、 向后差分格式程序function diffusion_2(h,lamda,time) %扩散方程向后差分格式，时间一阶精度，空间二阶精度，无条件稳定%h为空间步长，lamda为网格比，t为演化时间%h=0.1;time=1;lamda=0.5;tao=h2*lamda;t=0:tao:time;x=0:h:1;X,T=meshgrid(x,t);time1=time/tao+2; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点i=1:space1;j=1:time1;u(i,j)=0; %离散化的网格点wucha(1,i)=0;for i=1:space1-1 %边界条件初始化 u(i,1)=sin(pi*(i-1)*h);end%使用追赶法求每一时间层u的节点值,需要定义系数矩阵a和非其次项Zr=-lamda;m=1+2*lamda;kk=1:space1-2; %初始化系数矩阵 a(kk,kk)=0; a(1,1)=m; a(1,2)=r;for k=2:space1-3 %定义系数矩阵 a(k,k-1)=r; a(k,k)=m; a(k,k+1)=r;end a(space1-2,space1-3)=r; a(space1-2,space1-2)=m;%定义非齐次项 p=1:space1-2; q=1:time1-1; Z(p,q)=0; %矩阵a*x(i)=Z(i); %使用隐式格式计算 for jj=2:time1 %从第二时间层开始计算 %计算非齐次项 for m=1:space1-2 Z(m,jj-1)=u(m+1,jj-1); end %追赶法 E=chase(a,Z(:,jj-1); for kk=1:space1-2 u(kk+1,jj)=E(kk); endendu=u;figure(1);grid on;mesh(X,T,u);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(向后差分格式曲面图);%计算误差uu=exp(-pi2.*T).*sin(pi.*X);error=uu-u;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(向后差分格式误差曲面图);wucha=error(22,:)figure(3)plot(wucha,ro-);xlabel(t=0.1s，x轴);ylabel(误差轴);title(0.1秒时向后差分格式解与精确解的误差); 3、 Crank-Nicolson格式程序function diffusion_3(h,lamda,time) %扩散方程Crank-Nicolson差分格式，时间二阶精度，空间二阶精度，无条件稳定%h为空间步长，lamda为网格比，t为演化时间%h=0.1;time=1;lamda=0.5;tao=h2*lamda;t=0:tao:time;x=0:h:1;X,T=meshgrid(x,t);time1=time/tao+2; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点i=1:space1;j=1:time1;u(i,j)=0; %离散化的网格点wucha(1,i)=0;for i=1:space1-1 %边界条件初始化 u(i,1)=sin(pi*(i-1)*h);end%使用追赶法求每一时间层u的节点值,需要定义系数矩阵a和非其次项Zr=-0.5*lamda;r1=-r;m=1+lamda;m1=1-lamda;kk=1:space1-2; %初始化系数矩阵 a(kk,kk)=0; a(1,1)=m; a(1,2)=r;for k=2:space1-3 %定义系数矩阵 a(k,k-1)=r; a(k,k)=m; a(k,k+1)=r;end a(space1-2,space1-3)=r; a(space1-2,space1-2)=m;%定义非齐次项 p=1:space1-2; q=1:time1-1; Z(p,q)=0; %矩阵a*x(i)=Z(i); %使用隐式格式计算 for jj=2:time1 %从第二时间层开始计算 %计算非齐次项 for m=1:space1-2 Z(m,jj-1)=r1*u(m,jj-1)+m1*u(m+1,jj-1)+r1*u(m+2,jj-1); end %追赶法 E=chase(a,Z(:,jj-1); for kk=1:space1-2 u(kk+1,jj)=E(kk); end end u=u;figure(1);grid on;mesh(X,T,u);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(Crank-Nicolson差分格式曲面图);%计算误差uu=exp(-pi2.*T).*sin(pi.*X);error=uu-u;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel(x轴);ylabel(t轴);zlabel(u轴);title(Crank-Nicolson差分格式误差曲面图);wucha=error(22,:)figure(3)plot(wucha,ro-);xlabel(t=0.1s，x轴);ylabel(误差轴);title(0.1秒时Crank-Nicolson差分格式解与精确解的误差);4、 Du Fort-Frankel格式程序function diffusion_4(h,lamda,time) %扩散方程Du Fort-Frankel差分格式，时间二阶精度，空间二阶精度，无条件稳定%h为空间步长，lamda为网格比，t为演化时间%h=0.1;time=1;lamda=0.5;tao=h2*lamda;t=0:tao:time;x=0:h:1;X,T=meshgrid(x,t);time1=time/tao+2; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点i=1:space1;j=1:time1;u(i,j)=0; %离散化的网格点wucha(1,i)=0;for i=1:space1-1 %边界条