matlab解偏微分方程[1].pdf

《matlab解偏微分方程[1].pdf》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、第第第六六六章章章解解解偏偏偏微微微分分分方方方程程程1差差差分分分法法法解解解偏偏偏微微微分分分方方方程程程2热热热传传传导导导方方方程程程u=a2u的的的差差差分分分公公公式式式txx2.1显显显式式式格格格式式式热传导方程可以写成差分形式(右边取t时刻的值计算)u(x;t+¢t)¡u(x;t)u(x+¢x;t)¡2u(x;t)+u(x¡¢x;t)2¼a¢t(¢x)2¢t2即u(x;t+¢t)¼u(x;t)+a[u(x+¢x;t)¡2u(x;t)+u(x¡¢x;t)](¢x)2令x=i4x;t=j4t;i;j=0;1;2;¢¢¢n¡1;r=¢ta2,上式可写(¢x)2成显式差

2、分公式u(i;j+1)=(1¡2r)u(i;j)+r[u(i+1;j)¡2u(i;j)+u(i¡1;j)](¢x)2稳定条件为¢t·,截断误差为O((¢x)2;¢t).2a2例细杆传热问题定解问题是8>:u(x;t=0)='(x)其中0·x·20;a=10且(1;(10·x·11)'(x)=0;(x<10;x>11)根据上面的公式,编写如下程序x=0:20;t=0:0.01:1;a2=10;r=a2*0.01;u=zeros(21,101);u(10:11,1)=1;forj=1:100u(2:20,j+1)=(1-2*r)*

3、u(2:20,j)+r*(u(1:19,j)+u(3:21,j));plot(u(:,j));axis([02101]);pause(0.1)endsurf(u)2.2隐隐隐式式式格格格式式式热传导方程还可以写成如下差分形式(右边取t+4t时刻的值计算)u(x;t+¢t)¡u(x;t)u(x+¢x;t+¢t)¡2u(x;t+¢t)+u(x¡¢x;t+¢t)2¼a¢t(¢x)2引入与上面相同的足标,且令2ui+1;j+1¡2ui;j+1+ui¡1;j+1Hui;j+1=a(¢x)2则得到隐式公式为ui;j+1¡ui;j=Hui;j+1¢t令2ui+1;j¡2ui;j+ui¡1;jH

4、ui;j=a(¢x)2则显式公式写成ui;j+1¡ui;j=Hui;j¢t将显式与隐式相加,得平均公式ui;j+1¡ui;j11=Hui;j+Hui;j+1¢t22即是11+H¢t2ui;j+1=ui;j11¡H¢t2这个公式对任何步长都是恒稳定的.含时间的一维薛定谔方程一维运动的粒子的含时间的薛定格方程可以化成如下形式的抛物型问题(设方程右边的系数为1),@'@2'=i@t@x2利用前面推导的平均公式上面定义的算符H,得2311+iH4t627'i;j+1=45'i;j11¡iH4t2改写成下述形式23627'i;j+1=4¡15'i;j´Â¡'i;j11¡iH4t2中间函数Â

5、由下式定义，1(1¡iH4t)Â=2'i;j2即是(2i+H4t)Â=4i'i;j在使用MATLAB编程时,先计算出每一时刻的Ái;j对应的Â,程序中是把Â的系数表示为A,然后利用MATLAB的矩阵方程求解的功能求Â，再用它去求下一时刻的Ái;j+1.程序在一个220点的格子上，定义解析形式的位势(方形势阱或势垒，Gauss势阱或势垒，位势台阶或抛物线势阱)，建立一个初始的Gauss波包或Lorentz波包，它具有规定的平均位置、动量和空间宽度，并在保持格子的端点上'为零的边界条件下演化．用动画在每一时间步长上都显示概率密度j'j2和位势，以及波包在一指定点的左边和右边两部分每―

6、部分的总概率和平均位置．程序使用的数据是，选位势为一个方形位垒，高度是0.18,位于格点编号为105和115的位置。初始的波包是一高斯波包，其平均波数和宽度为k0=0:6;¾=10,中心位置x0在编号为40的格点上。时间步长为1，演化时间为130。动画演示了波包向势垒行进，在势垒上发生透射和反射，最后分为透射波和反射波二部分向不同的方向运动。下图给出了这个过程中的几个画面。NPTS=220;sigma=10;k0=0.6;x0=40;time=130;v(NPTS)=0;v(105:115)=+0.18;A=diag(-2+2i+v)+diag(ones(NPTS-1,1),1)

7、+diag(ones(NPTS-1,1),-1);PHI0(1)=0;PHI0(NPTS)=0;forx=2:NPTS-1;PHI0(x)=exp(0.5i*x)*exp(-(x-x0)^2*log10(2)/sigma^2);endPHI(:,1)=PHI0.';fory=2:timeCHI(:,y)=4i*(APHI(:,y-1));PHI(:,y)=CHI(:,y)-PHI(:,y-1);endmo=moviein(time);forj=1:timeplot([1:NPT