维波动方程的有限差分法

1、维波动方程的有限差分法-作者 : _-日期 : _学 生 实 验 报 告实验课程名称偏微分方程数值解开课实验室数统学院学 院数 统 年级 2013 专业班 信计 02 班学 生 姓 名学 号开 课 时 间 2015至 2016学年第2 学期总 成 绩教师签名2数学与统计学院制开课学院、实验室： 数统学院实验时间 ：2016年 6 月 20 日实验项目实验项目类型一维波动方程的有限差分法名 称验证演示综合设计其他指导教成 绩师曾芳是一实验目的通过该实验，要求学生掌握求解一维波动方程的有限差分法，并能通过计算机语言编程实现。二实验内容考虑如下的初值问题：2u2u, x0,1 , t0, 2t2x2

2、ux,0sinx, ux,00 (1)tu0, tu1,t0, t0,21在第三部分写出问题（1）三层显格式。2根据你写出的差分格式，编写有限差分法程序。将所写程序放到第四部分。3取 h0.1,0.1h ，分别将 t0.5,1.0,1.5,2.0 时刻的数值解画图显示。4. 该问题的解析解为 u x,t cos t sin x ，将四个时刻的数值解的误差画图显示，对数值结果进行简单的讨论。三实验原理、方法（算法）、步骤1、三层显格式建立由于题中 h0.1,0.1h ， x0,1 , t0,2 ，取 n10, m 200 ，故令网比r0.1 ， xjj h, j0,1,2,l10， tkk,k0

3、,1,l, 200 ，在 x j , tk 内网个点处，利用二h阶中心差商得到如下格式：ukj12ukjukj1o2ukj 12ukju kj 1o h2（ 2）2h2略去误差项得到：3ukj1r 2ukj 12 1 r 2 ukjr 2 ukj 1ukj1（3）其中 j1,2,l 9, k1,2,l,199 ，局部截断误差为 o2h2。对于初始条件 ux,0sinx ，建立差分格式为：u0jsinxjsinjh , j0,1,2,l 10（ 4）对于初始条件tu x,00 ，利用中心差商，建立差分格式为：u1j1uj0, 即 u1j =u j1， j0,1,2,l10（ 5）2对于边界条件

4、u0,tu 1, t0, t0,2 ，建立差分格式为：u0kunk0, k0,1,l, 200（ 6）将差分格式延拓使 k0 为内点，代入（ 3）得到的式子再与（5）联立消去 u j1后整理得到：u1j1 r 2u0j 11 r 2 u 0j1 r 2 u0j 1（ 7）22综上（ 3）、（ 4）、（ 6）、（ 7）得到三层显格式如下：（局部截断误差为 o2h2 ）k 1r2 k2 1r2k2 kk 1, jl9, k1,2,l,199u ju j 1u jr uj1u j1,2,0sinx jsinjh, j0,1,2,lu j10（8）u1j1 r 2u0j 11r 2u 0j1 r 2u

5、0j1, j 1,2,l92u0kunk20,k0,1,l, 200其中 r0.1。h四实验环境（所用软件、硬件等）及实验数据文件matlab三层显格式程序如下：%一维波动方程，三层显格式求解法h=0.1;tau=0.1*h;r=tau/h;n=1/h;m=2/tau;x=0:h:1;t=0:tau:2;u=sin(pi*x);%计算 t=0 时刻的 u 值u(1,11)=0;for j=2:nu(2,j)=0.5*r2*u(1,j+1)+(1-r2)*u(1,j)+0.5*r2*u(1,j-1);end%定义 x=0 边界上的数值for k=1:m+1u(k,1)=0;end4%定义 x=1

6、 边界上的数值for k=1:m+1u(k,n+1)=0;end%迭代计算开始 , 差分格式for k=2:mfor j=2:nu(k+1,j)=r2*u(k,j+1)+2*(1-r2)*u(k,j)+r2*u(k,j-1)-u(k-1,j);endendu(201,:)=zeros(1,11);%计算 k=201 行的数值解u2(201,11)=0;for j=2:nu2(201,j)=r2*u(200,j+1)+2*(1-r2)*u(200,j)+r2*u(200,j-1)-u(199,j);endu=u+u2;u=rot90(u,2);%将矩阵 u 旋转 180 度赋值于 u%作出图像x

7、,t=meshgrid(0:0.1:1,0:0.01:2);%划分网格%作出数值解的函数图像subplot(2,2,1);mesh(x,t,u);title(u(x,t)数值解的函数图像 );xlabel(x变量 );ylabel(t变量 );zlabel(u值);%作出精确解的函数图像subplot(2,2,2);u1=cos(pi*t).*sin(pi*x);mesh(x,t,u1);title(u(x,t)精确解的函数图像 );xlabel(x变量 );ylabel(t变量 );zlabel(u值);%作出 t=0.5 ，1.0 ，1.5, 2.0时刻的绝对误差图像5subplot(2,

8、2,3);wucha=abs(u-u1);x=0:h:1;plot(x,wucha(51,:),g*-);hold ongrid onplot(x,wucha(101,:),ro-);hold onplot(x,wucha(151,:),ks-);hold onplot(x,wucha(201,:),mp-);title(t=0.5， 1.0 ，1.5, 2.0时刻的绝对误差函数图像);xlabel(x变量 ); ylabel(绝对误差值 );legend(t=0.5,t=1.0,t=1.5,t=2.0);%作出 t=0.5 ，1.0 ，1.5, 2.0时刻的数值解函数图像subplot(2,

9、2,4);x=0:h:1;plot(x,u(51,:),g*-);hold ongrid onplot(x,u(101,:),ro-);hold onplot(x,u(151,:),ks-);hold onplot(x,u(201,:),mp-);title(t=0.5， 1.0 ，1.5, 2.0时刻的数值解函数图像 );xlabel(x变量 ); ylabel(u值 );legend(t=0.5,t=1.0,t=1.5,t=2.0);%当然也可以作出u(x,t)绝对误差的函数图像%mesh(x,t,wucha);%title(u(x,t)绝对误差的函数图像 );%xlabel(x变量 );

10、%ylabel(t变量 );%zlabel(绝对误差值 );6五实验结果及实例分析1、 u(x,t)在 t=0.5,1.0,1.5,2.0时刻的数值解、精确解以及绝对误差表 1 u(x,t)在 t=0.5,1.0,1.5,2.0时刻的数值解时刻 tt=0.5,1.0,1.5,2.0时刻的数值解t=0.50 -0.0059-0.0113-0.0155-0.0182-0.0192-0.0182-0.0155-0.0113-0.00590t=1.00 -0.3090-0.5877-0.8090-0.9510-0.9999-0.9510-0.8090-0.5877-0.30900t=1.500.002

11、00.00380.00520.00610.00640.00610.00520.00380.00200t=2.000.30900.58780.80900.95111.00000.95110.80900.58780.30900表 2 u(x,t)在 t=0.5,1.0,1.5,2.0时刻的精确解时刻 tt=0.5,1.0,1.5,2.0时刻的精确解t=0.500.00000.00000.00000.00000.00000.00000.00000.00000.00000t=1.00-0.3090 -0.5878 -0.8090-0.9511 -1.0000 -0.9511 -0.8090 -0.58

12、78 -0.3090 0t=1.500.00000.00000.00000.00000.00000.00000.00000.00000.00000t=2.000.30900.58780.80900.95111.00000.95110.80900.58780.30900表 3 u(x,t)在 t=0.5,1.0,1.5,2.0时刻的绝对误差时刻 tt=0.5,1.0,1.5,2.0时刻的绝对误差t=0.500.00590.01130.01550.01820.01920.01820.01550.01130.00590t=1.000.00000.00000.00010.00010.00010.00010.00010.00000.00000t=1.500.00200.00380.00520.00610.00640.00610.00520.00380.00200t=2.000.00000.00000.00000.00000.00000.00000.00000.00000.00000说明：在 t=0.5 时刻的绝对误差最大，t=1.5 时刻次之， t=1 与 t=2 时刻的绝对误差均较小，由于 r0.11 ，该格式