《计算物理(本科)》[第10章]优质PPT.ppt

《计算物理(本科)》[第10章]优质PPT.ppt

《《计算物理(本科)》[第10章]优质PPT.ppt》由会员分享，可在线阅读，更多相关《《计算物理(本科)》[第10章]优质PPT.ppt（17页珍藏版）》请在冰豆网上搜索。

若弦线是均匀，即迫力。

若弦线是均匀，即（x）=为为常数，则称波动方程为常数，则称波动方程为一维非齐次波动方程，即一维非齐次波动方程，即二维波动方程为二维波动方程为三维波动方程为三维波动方程为其中其中v=（T/）1/2称波速。

称波速。

合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法10.2一维波动方程的差分法一维波动方程的差分法一、波动方程：

一、波动方程：

均匀弦线的自由振动的波动方程为均匀弦线的自由振动的波动方程为其中其中v波速；

波速；

l为弦线长度；

为弦线长度；

T为有限时间。

为有限时间。

t=tk=kk=0,1,2,M;

=T/M（时间步序号时间步序号）x=xi=ihi=0,1,N;

h=l/N（空间步序号）空间步序号）建立差分格式建立差分格式同上章求热传导方程，构造网同上章求热传导方程，构造网格节点，如图所示格节点，如图所示txk+1ki+1ii-1MN合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法用二用二阶阶中心差商作为二阶微商近似得中心差商作为二阶微商近似得由上两式，并由上两式，并令令=v/h可得可得一维波动方程的差分格式一维波动方程的差分格式其中其中N=l/h,M=T/,“”表示取整。

表示取整。

合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法几何关系几何关系各点关系如图所示，从图中直观各点关系如图所示，从图中直观地看出，地看出，k+1时刻第时刻第i个点的值是个点的值是由由k时刻的第时刻的第i-1、i、i+1三点和三点和k-1时刻的时刻的i的值推算出来的。

的值推算出来的。

txk+1ki+1ii-1MN二、边界条件二、边界条件将其将其变为差分格式为变为差分格式为合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法三、初始条件三、初始条件（一般形式）（一般形式）上式使用中心一阶差商近似的差分格式为上式使用中心一阶差商近似的差分格式为合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法再再利用上面一维波动方程的差分格式利用上面一维波动方程的差分格式令令k=0得：

得：

结合结合得：

合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法四、归纳上述一维波动方程的差分格式为四、归纳上述一维波动方程的差分格式为上式的差分上式的差分格式收敛稳定条件为格式收敛稳定条件为=v/h1合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法10.3MATLAB实现波动方程数值解法实现波动方程数值解法MATLAB提供了一个求解波动方程的提供了一个求解波动方程的hyperbolic（）函数，函数，此函数主要是用有限元算法进行编写程序的。

此函数主要是用有限元算法进行编写程序的。

其用法为其用法为u1=hyperbolic（u0,ut0,tlist,b,p,e,t,c,a,f,d）其中其中u0为为初始值；

初始值；

ut0为为u一阶导数初始值；

一阶导数初始值；

tlist为时间序为时间序列；

列；

b是边界条件，可以用矩阵形式，也可以用是边界条件，可以用矩阵形式，也可以用m文件格文件格式，主要依赖时间式，主要依赖时间t；

p,e,t为网格参数；

为网格参数；

c,a,f,d是方程系是方程系数，可以是时间的函数；

另外可以在函数末加入误差限数，可以是时间的函数；

另外可以在函数末加入误差限来控制。

对于标量形式的波动方程，函数返回值为一个来控制。

对于标量形式的波动方程，函数返回值为一个矩阵。

矩阵。

设波动方程形式为设波动方程形式为合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法解解：

在：

在MATLAB中建立中建立m文件文件jswlx_10_3_1.m【例例1】求波动方程求波动方程求解的求解的范围为正方形区域：

范围为正方形区域：

-1x,y1，t0；

初始条件：

当初始条件：

当-1x,y1时时,第一类第一类u（x,y,0）=atancos（2x）第二类第二类du（x,y,0）/dt=3sin（x）expsin（y/2）；

边界条件：

当边界条件：

当t0时，时，u（-1,y,t）=0，u（1,y,t）=0。

合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法%求解求解标准波动方程标准波动方程d2u/dt2-div（grad（u）=0%g=squareg;

%定义单位方形区域函数定义单位方形区域函数b=squareb3;

%定义左右定义左右0边界条件边界条件,顶底顶底0导数边界条件函数导数边界条件函数c=1;

a=0;

f=0;

d=1;

p,e,t=initmesh（squareg）;

%网格化网格化%定义初始条件，选择第一类初始条件定义初始条件，选择第一类初始条件u（0）=atan（cos（pi/2*x）和和%第二类初始条件第二类初始条件dudt（0）=3*sin（pi*x）.*exp（sin（pi/2*y）,以免把大量的精力以免把大量的精力%花费在高阶振动模型上花费在高阶振动模型上,这样可以得到可行的时间步长这样可以得到可行的时间步长.实现如下实现如下:

x=p（1,:

）;

y=p（2,:

u0=atan（cos（pi/2*x）;

ut0=3*sin（pi*x）.*exp（sin（pi/2*y）;

n=31;

tlist=linspace（0,5,n）;

%在时间段为在时间段为0,5的的31个点上求解个点上求解uu=hyperbolic（u0,ut0,tlist,b,p,e,t,c,a,f,d）;

%求此波动求此波动（双曲型双曲型）方程问题方程问题%矩形网格插值矩形网格插值.为了加速绘图为了加速绘图,在矩形网格下插值在矩形网格下插值delta=-1:

0.1:

1;

uxy,tn,a2,a3=tri2grid（p,t,uu（:

1）,delta,delta）;

gp=tn;

a2;

a3;

合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法%动画图形显示动画图形显示newplot;

M=moviein（n）;

umax=max（max（uu）;

umin=min（min（uu）;

fori=1:

n,.ifrem（i,10）=0,.end,.pdeplot（p,e,t,xydata,uu（:

i）,zdata,uu（:

i）,zstyle,continuous,.mesh,off,xygrid,on,gridparam,gp,colorbar,off）;

.axis（-11-11uminumax）;

caxis（uminumax）;

.M（:

i）=getframe;

.ifi=n,.fprintf（donen）;

.end,.endnfps=5;

movie（M,2,nfps）;

%演示动画演示动画合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法解解：

在MATLAB中建立中建立m文件文件jswlx_10_3_2.m求解的求解的范围为正方形区域：

当-1x,y1时时,第一类第一类u（x,y,0）=sin（1+x）/2第二类第二类du（x,y,0）/dt=（1+x）（1-x）；

【例例2】

（课本课本P117）求波动方程求波动方程合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法%求解求解标准波动方程标准波动方程d2u/dt2-4div（grad（u）=0%g=squareg;

%网格化网格化%定义初始条件，选择第一类初始条件定义初始条件，选择第一类初始条件u0=sin（pi/2*（1+x）和和%第二类初始条件第二类初始条件dudt（0）=1/4*（1+x）.*（1-x）,以免把大量的精力以免把大量的精力%花费在高阶振动模型上花费在高阶振动模型上,这样可以得到可行的时间步长这样可以得到可行的时间步长.实现如下实现如下:

u0=sin（pi/2*（1+x）;

ut0=1/4*（1+x）.*（1-x）;

%演示动画演示动画合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法实验九实验九波动方程和薛定谔方程求解波动方程和薛定谔方程求解一、实验目的一、实验目的1了解一维波动方程的差分解法；

了解一维波动方程的差分解法；

2掌握二维波动方程的差分解法。

掌握二维波动方程的差分解法。

二、实验内容二、实验内容1建立一维波动方程的差分格式；

建立一维波动方程的差分格式；

2建立二维波动方程的差分格式；

建立二维波动方程的差分格式；

2会用会用MATLAB求解波动方程的求解波动方程的hyperbolic（）函数求解函数求解某些波动方程。

某些波动方程。

合肥工业大学电子科学与应用物理学院第第十十章章波波动动方方程程的的数数值值解解法法三、练习题三、练习题1.在在Matlab中，用求解波动方程函数中，用求解波动方程函数u1=hyperbolic（u0,ut0,tlist,b,p,e,t,c,a,f,d）求解下面问题，并绘图展示变化过程：

求