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

一维波动方程的有限差分法学 生 实 验 报 告实验课程名称 偏微分方程数值解 开课实验室 数统学院 学 院 数 统 年级 2013 专业班 信计02班 学 生 姓 名 学 号 开 课 时 间 2015 至 2016 学年第 2 学期总 成 绩教师签名数学与统计学院制开课学院、实验室： 数统学院 实验时间 ： 2016年 6月20日实验项目名 称一维波动方程的有限差分法实验项目类型验证演示综合设计其他指导教师曾芳成 绩是一实验目的通过该实验，要求学生掌握求解一维波动方程的有限差分法，并能通过计算机语言编程实现。二实验内容考虑如下的初值问题： (1)1在第三部分写出问题（1）三层显格式。2根据你写出的差分格式，编写有限差分法程序。将所写程序放到第四部分。3取，分别将时刻的数值解画图显示。4. 该问题的解析解为，将四个时刻的数值解的误差画图显示，对数值结果进行简单的讨论。三实验原理、方法（算法）、步骤1、三层显格式建立由于题中，取，故令网比，在内网个点处，利用二阶中心差商得到如下格式： （2）略去误差项得到： （3）其中，局部截断误差为。对于初始条件，建立差分格式为： （4）对于初始条件，利用中心差商，建立差分格式为： （5）对于边界条件，建立差分格式为： （6）将差分格式延拓使为内点，代入（3）得到的式子再与（5）联立消去后整理得到： （7）综上（3）、（4）、（6）、（7）得到三层显格式如下：（局部截断误差为） （8）其中。四实验环境（所用软件、硬件等）及实验数据文件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:N u(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+1 u(k,1)=0;end%定义x=1边界上的数值for k=1:M+1 u(k,N+1)=0;end%迭代计算开始,差分格式for k=2:M for j=2:N u(k+1,j)=r2*u(k,j+1)+2*(1-r2)*u(k,j)+r2*u(k,j-1)-u(k-1,j); endend u(201,:)=zeros(1,11);%计算k=201行的数值解u2(201,11)=0;for j=2:N u2(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,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时刻的绝对误差图像subplot(2,2,3);wucha=abs(u-u1);x=0:h:1;plot(x,wucha(51,:),g*-);hold on grid 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,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变量);%ylabel(t变量);%zlabel(绝对误差值);五实验结果及实例分析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.0059 0 t=1.00 -0.3090 -0.5877 -0.8090 -0.9510 -0.9999 -0.9510 -0.8090 -0.5877 -0.3090 0 t=1.50 0.0020 0.0038 0.0052 0.0061 0.0064 0.0061 0.0052 0.0038 0.0020 0 t=2.00 0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0 表2 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.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 t=1.00 -0.3090 -0.5878 -0.8090 -0.9511 -1.0000 -0.9511 -0.8090 -0.5878 -0.3090 0 t=1.50 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 t=2.00 0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0 表3 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.0059 0 t=1.00 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0001 0.0000 0.0000 0 t=1.50 0.0020 0.0038 0.0052 0.0061 0.0064 0.0061 0.0052 0.0038 0.0020 0 t=2.00 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 说明：在t=0.5时刻的绝对误差最大，t=1.5时刻次之，t=1与t=2时刻的绝对误差均较小，由于，该格式稳定，由数值计算得到的矩阵不难看出，数值解符合理论解。2、u(x,t)在t=0.5,1.0,1.5,2.0时刻的数值解、绝对误差函数图像 图1 数值解、精确解以及绝对误差函数图像说明：上两图为函数的数值解与精确解，下两图为 t=0.5,1.0,1.5,2.0时刻的数值解、绝对误差函数图像，符合理论解。教师签名年 月 日11