一维波动方程的有限差分法

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

2、问题（1）三层显格式。2.根据你写出的差分格式，编写有限差分法程序。将所写程序放到第四部分。3.取h 0.1,0.1h，分别将t 0.5,1.0,1.5,2.0时刻的数值解画图显示。4该问题的解析解为u x,t cos tsin x，将四个时刻的数值解的误差画图显示,对数值结果进行简单的讨论。三.实验原理、方法（算法）、步骤1、三层显格式建立由于题中h 0.1,0.1h ,Xj j h, j 0,1,2,L 10 , tkkx 0,1 ,t 0,2，取 N 10, M 200，故令网比 r 0.1 ,h,k 0,1,L ,200，在Xj,tk内网个点处，利用二阶中心差商得到如下格式:对于边界条

3、件u 0,t u 1,t0,t0,2，建立差分格式为：u； uN 0,k 0,1,L ,200（6）将差分格式延拓使k 0为内点，代入（3）得到的式子再与（5）联立消去u1后整理得到:1Uj2印11r2 uq 1j 22 0r Uj 1(7)综上（3）、（4）、(6)、(7)得到三层显格式如下：（局部截断误差为o 2 h2 ）k 1Uj2 kr Uj 1 22k1 r Uj2 kk 1r Uj 1 Uj ,j1,2,L 9,k 1,2,L ,1990Ujsinxjsin jh , j0,1,2,L 1011Ujj22 0r Uj 112 0 1 2 0r Uj - r Uj 1, j 1,2,

4、L 9(8)kkUqUn0,k0,1L,200略去误差项得到:其中j 1,2丄9,k对于初始条件对于初始条件k 1k k 1Uj2u j Ujkk kUj 1 2Uj Uj 1h2o h2k 1UjkUj2 kr Uj2 ox，建立差分格式为：sin xjsin jh , j2 112 L ,199，局部截断误差为 u x,0 sin0Ujk 1Ujh2 。0,1,2,L 10一 u x,0 t利用中心差商，建立差分格式为:1Uj21Uj0,即u1二Uj1, j 0,1,2丄 10(2)(3)(4)(5)其中r 0.1 h四实验环境（所用软件、硬件等）及实验数据文件Matlab三层显格式程序如

5、下:%维波动方程，三层显格式求解法 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*rA2*u(1,j+1)+(1-rA2)*u(1,j)+0.5*rA2*u(1,j-1);end淀义x=0边界上的数值for k=1:M+1u(k,1)=0;end淀义x=1边界上的数值for k=1:M+1u(k,N+1)=0;end%迭代计算开始，差分格式for k=2:Mfor j=2:N u(k+1,j)=rA2*u(k,j+1)+2

6、*(1-rA2)*u(k,j)+rA2*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)=rA2*u(200,j+1)+2*(1-rA2)*u(200,j)+rA2*u(200,j-1)-u(199,j);endu=u+u2;u=rot90(u,2);% 将矩阵u旋转180度赋值于u%乍出图像x,t=meshgrid(0:0.1:1,0:0.01:2);%划分网格%乍出数值解的函数图像subplot(2,2,1);mesh(x,t,u);title(u(x,t)数值

7、解的函数图像);xlabel(x 变量);ylabel(t 变量);zlabel(u 值);%乍出精确解的函数图像subplot(2,2,2);u1=cos(pi*t).*s in (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时刻的绝对误差图像subplot(2,2,3);wucha=abs(u-u1);x=0:h:1;plot(x,wucha(51,:),g*-);hold ongrid onplot(x,wucha(101,:

8、),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,2,4);x=0:h:1;plot(x,u(51,:),g*-);hold ongrid onplot(x,u(101,:),ro-);hold onplot(x,u

9、(151,:),ks-);hold onplot(x,u(201,:),mp-);title(t=0.5 , 1.0 , 1.5, 2.0时刻的数值解函数图像);xlabel(x 变量);ylabel(u值);lege nd(t=0.5,t=1.0,t=1.5,t=2.0);%i然也可以作出u(x,t)绝对误差的函数图像%mesh(x,t,wucha);%title(u(x,t)绝对误差的函数图像);%xlabel(x 变量);%ylabel(t 变量);%zlabel(绝对误差值);1、u(x,t) 在 t=0.5,1.0,1.5,2.0表 1 u(x,t)五实验结果及实例分析时刻的数值解、

10、精确解以及绝对误差在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.0059 0t=1.00-0.309005877 -0.8090-0.9510-0.9999-0.9510 -0.8090-0.5877-0.3090 0t=1.500.00200.00380.00520.00610.00640.00610.00520.00380.00200t=2.000.30900.58780.80900.95111

11、.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.0000 0.00000.00000.00000.00000.00000t=1.00-0.3090 -0.5878 -0.8090-0.9511-1.0000 -0.9511 -0.8090 -0.5878 -0.3090 0t=1.500.00000.00000.00000.0000 0.00000.00000.00000.00000.00000t=2.000.3

12、0900.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.0001 0.00010.00010.00010.00000.00000t=1.500.00200.00380.00520.00610.00640.00610.00520.00380.00200t=2.000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000说明：在t=0.5时刻的绝对误差最大，t=1.5时