偏微分方程程序附件Word文档格式.docx

1、for n=1:5/tao-1 for m=1:10/h-1 U（m,n）=u（m+1,n+1）;U=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: for mm=1: uu（mm,nn）=（nn*tao）2*sin（mm*h）;UU=uuerror= U-UU;figure（2）;mesh（X,T,error）;%mesh（X,T,UU）;迎风格式误差曲面图figure（3）;mesh（

2、X,T,UU）;古典解2、 Lax-Wendroff格式function lax_w1（h,lamda） %Lax-Wendroff 二阶精度 显示格式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

3、）+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差分格式Lax-Wendroff格式差分曲面图error=U-UU;Lax-Wendroff格式误差曲面图3、 隐式中心格式function implicit_1（h,lamda） %隐式格式，精度为：时间一阶，空间二阶；无条件稳定time1=5/tao; %时间为1时的时间层网格点space1=10/h; %空间为1时的时间

4、层网格点space1;time1;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; a（space1-1,space1-1）=1; p=1: q=1:time1-1; Z（p,q）=0; %矩阵a*x（i）=Z（i）;time1space1-1 Z（m-1,n-1）=u（m,n-1）+tao*（2*（n-1）*tao*sin（m-1）*h）+. n2*tao2*cos（m

5、*h）; 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: u（v+1,n）=E（v）;隐式格式差分曲面图error=UU-U;隐式格式误差曲面图调用追赶法接三对角矩阵函数chase.m为：function X=chase（A,d）%追赶法只用来解三对角方程组，其中A为方程

6、的系数，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）;X（n）=y（n）/U（n）;for i=n-1:-1:1 X（i）=（y（i）-c（i）*X（i+1）/U（i）;二、 计算下面双曲方程的近似解1、 二阶显示格式function explicit_2（h,lamda

7、） %波动方程显示格式，二阶精度，lamda1时稳定tao=h*lamda;1-tao;1-h;time1=1/tao+1;space1=1/h+1;time1 %边界条件离散 u（space1,n）=sin（n-1）*tao）;for 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）;for n=3:1/taospace1-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）+. tao

8、2*（n-1）*tao）2-（m-1）*h）2）*sin（n-1）*tao*（m-1）*h）; %波动方程显示格式time1-2space1-2波动方程显示格式差分曲面图%计算古典解space1-2;time1-2; uu（mm,nn）=sin（mm*h*nn*tao）;波动方程显示格式误差曲面图波动方程古典解2、 隐式二阶格式function implicit_2（h,lamda） %波动方程隐式格式，二阶精度，无条件稳定 %离散化的网格点v（i,j）=0; %定义关于时间的一阶导数for i=1: v（i,1）=（i-1）*h; %时间的一阶导数初始值离散化 for j=1: w（i,j）

9、=0; %定义关于空间的一阶导数end w（i,1）=0; %空间的一阶导数初始值离散化%使用追赶法求每一时间层v的节点值,需要定义系数矩阵a和非其次项Zr=-0.25*lamda2;m=1+0.5*lamda2; a（1,1）=m;space1-3 %定义系数矩阵 a（k,k-1）=r; a（k,k）=m; a（space1-2,space1-3）=r; a（space1-2,space1-2）=m;%定义非齐次项 %使用隐式格式计算 for jj=2:time1 %从第二时间层开始计算 v（1,jj）=0; %v左边界值离散 v（space1,jj）=（u（space1,jj）-u（spa

10、ce1,jj-1）/tao; %v右边界值离散 %计算非齐次项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）; 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）;,jj-1）; for kk=1:

11、v（kk+1,jj）=E（kk）; %计算w（:,jj）和u（:,jj） for ii=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）; for g=2: u（g,jj）=u（g,jj-1）+tao*v（g,jj）;i1=1:uu（i1,jj）=0;三、 计算下面抛物方程的近似解1、 向前差分格式程序：function diffusion_1（h,lamda,time） %扩散方程向前差分格式，时间一阶精度，空间二阶精度，lamda=0,5稳定%h为空间步长，l

12、amda为网格比，t为演化时间%h=0.1;time=1;lamda=0.5;tao=h2*lamda;t=0:time;x=0:1;time1=time/tao+2;wucha（1,i）=0;space1-1 %边界条件初始化 u（i,1）=sin（pi*（i-1）*h）; u（m,n）=（1-2*lamda）*u（m,n-1）+lamda*u（m-1,n-1）+lamda*u（m+1,n-1）;u=umesh（X,T,u）;向前差分格式曲面图uu=exp（-pi2.*T）.*sin（pi.*X）;error=uu-u;向前差分格式误差曲面图%plot（error）;mesh（X,T,uu）

13、;wucha=error（22,:）figure（4）plot（wucha,ro-t=0.1s，x轴误差轴0.1秒时差分格式解与精确解的误差2、 向后差分格式程序function diffusion_2（h,lamda,time） %扩散方程向后差分格式，时间一阶精度，空间二阶精度，无条件稳定%使用追赶法求每一时间层u的节点值,需要定义系数矩阵a和非其次项Zr=-lamda;m=1+2*lamda;kk=1:a（kk,kk）=0;time1 %从第二时间层开始计算 Z（m,jj-1）=u（m+1,jj-1）; %追赶法 u（kk+1,jj）=E（kk）;向后差分格式曲面图向后差分格式误差曲面图figure（3）0.1秒时向后差分格式解与精确解的误差3、 Crank-Nicolson格式程序function diffusion_3（h,lamda,time） %扩散方程Crank-Nicolson差分格式，时间二阶精度，空间二阶精度，无条件稳定r=-0.5*lamda;r1=-r;m=1+lamda;m1=1-lamda;%使用隐式格式计算 for jj=2: Z（m,jj-1）=r1*u（m,jj-1）+m1*u（m+1,jj-1）+r1*u（m+2,jj-1）; end