资源描述:
《一维波动方程的有限差分法》由会员上传分享,免费在线阅读,更多相关内容在应用文档-天天文库。
1、学生实验报告实验课程名称偏微分方程数值解开课实验室数统学院学院数统年级2013专业班信计02班学生姓名学号开课时间2015至2016学年第2学期总成绩教师签名数学与统计学院制开课学院、实验室:数统学院实验时间:2016年6月20日实验项目名称一维波动方程的有限差分法实验项目类型验证演示综合设计其他指导教师曾芳成绩是一.实验目的通过该实验,要求学生掌握求解一维波动方程的有限差分法,并能通过计算机语言编程实现。二.实验内容考虑如下的初值问题:(1)1.在第三部分写出问题(1)三层显格式。2.根据你写出的差分格式,编写有限差分法程序。将所写程序放到第四部分。3.取,分别将时刻的数
2、值解画图显示。4.该问题的解析解为,将四个时刻的数值解的误差画图显示,对数值结果进行简单的讨论。三.实验原理、方法(算法)、步骤1、三层显格式建立由于题中,,取,故令网比,,,在内网个点处,利用二阶中心差商得到如下格式:(2)略去误差项得到:(3)其中,局部截断误差为。对于初始条件,建立差分格式为:(4)对于初始条件,利用中心差商,建立差分格式为:(5)对于边界条件,建立差分格式为:(6)将差分格式延拓使为内点,代入(3)得到的式子再与(5)联立消去后整理得到:(7)综上(3)、(4)、(6)、(7)得到三层显格式如下:(局部截断误差为)(8)其中。四.实验环境(所用软件、
3、硬件等)及实验数据文件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;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;end%定义x=1边界上的数值fork=1:M+1u(k,N+1)=0;end%迭代计算开始,差分格式fork=2:Mf
4、orj=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,:)=zeros(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);%划分网格%作出数值解的函数图像s
5、ubplot(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;plo
6、t(x,wucha(51,:),'g*-');holdongridonplot(x,wucha(101,:),'ro-');holdonplot(x,wucha(151,:),'ks-');holdonplot(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;p
7、lot(x,u(51,:),'g*-');holdongridonplot(x,u(101,:),'ro-');holdonplot(x,u(151,:),'ks-');holdonplot(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);%t