资源描述：
《维波动方程的有限差分法.docx》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、维波动方程的有限差分法--------------------------------------------------------------------------作者:_____________--------------------------------------------------------------------------日期:_____________学生实验报告实验课程名称偏微分方程数值解开课实验室数统学院学院数统年级2013专业班信计02班学生姓名学号开课时间2015至2016学年第2学期总成绩教师签名2数学
2、与统计学院制开课学院、实验室：数统学院实验时间：2016年6月20日实验项目实验项目类型一维波动方程的有限差分法名称验证演示综合设计其他指导教成绩师曾芳是一．实验目的通过该实验，要求学生掌握求解一维波动方程的有限差分法，并能通过计算机语言编程实现。二．实验内容考虑如下的初值问题：2u2u,x0,1,t0,2t2x2ux,0sinx,ux,00(1)tu0,tu1,t0,t0,21．在第三部分写出问题（1）三层显格式。2．根据你写出的差分格式，编写有限差分法程序。将所写程序放到第四部分。3．取h0.1,0.1h，分别将t0.5,1.0,1.5
3、,2.0时刻的数值解画图显示。4.该问题的解析解为ux,tcostsinx，将四个时刻的数值解的误差画图显示，对数值结果进行简单的讨论。三．实验原理、方法（算法）、步骤1、三层显格式建立由于题中h0.1,0.1h，x0,1,t0,2，取N10,M200，故令网比r0.1，xjjh,j0,1,2,L10，tkk,k0,1,L,200，在xj,tk内网个点处，利用二h阶中心差商得到如下格式：ukj12ukjukj1o2ukj12ukjukj1oh2（2）2h2略去误差项得到：3ukj1r2ukj121r2ukjr2ukj1ukj1（3）其中j1
4、,2,L9,k1,2,L,199，局部截断误差为o2h2。对于初始条件ux,0sinx，建立差分格式为：u0jsinxjsinjh,j0,1,2,L10（4）对于初始条件tux,00，利用中心差商，建立差分格式为：u1j1uj0,即u1j=uj1，j0,1,2,L10（5）2对于边界条件u0,tu1,t0,t0,2，建立差分格式为：u0kuNk0,k0,1,L,200（6）将差分格式延拓使k0为内点，代入（3）得到的式子再与（5）联立消去uj1后整理得到：u1j1r2u0j11r2u0j1r2u0j1（7）22综上（3）、（4）、（6）、（
5、7）得到三层显格式如下：（局部截断误差为o2h2）k1r2k21r2k2kk1,jL9,k1,2,L,199ujuj1ujruj1uj1,2,0sinxjsinjh,j0,1,2,Luj10（8）u1j1r2u0j11r2u0j1r2u0j1,j1,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
6、(pi*x);%计算t=0时刻的u值u(1,11)=0;forj=2:Nu(2,j)=0.5*r^2*u(1,j+1)+(1-r^2)*u(1,j)+0.5*r^2*u(1,j-1);end%定义x=0边界上的数值fork=1:M+1u(k,1)=0;end4%定义x=1边界上的数值fork=1:M+1u(k,N+1)=0;end%迭代计算开始,差分格式fork=2:Mforj=2:Nu(k+1,j)=r^2*u(k,j+1)+2*(1-r^2)*u(k,j)+r^2*u(k,j-1)-u(k-1,j);endendu(201,:)=zer
7、os(1,11);%计算k=201行的数值解u2(201,11)=0;forj=2:Nu2(201,j)=r^2*u(200,j+1)+2*(1-r^2)*u(200,j)+r^2*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变量');ylabe
8、l('t变量');zlabel('u值');%作出精确解的函数图像subplot(2,2,2);u1=cos(pi*t).*sin(pi*x);mesh(x,t,u1);titl
