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

学生实验报告实验课程名称偏微分方程数值解开课实验室数统学院学院数统年级2013专业班信计02班学生姓名学号开课时间2015至2016学年第2学期总成绩教师签名数学与统计学院制2开课学院、实验室：数统学院实验时间：2016年6月20日实验项目名称一维波动方程的有限差分法实验项目类型验证演示综合设计其他指导教师曾芳成绩是一．实验目的通过该实验，要求学生掌握求解一维波动方程的有限差分法，并能通过计算机语言编程实现。二．实验内容考虑如下的初值问题：2222,0,1,0,2,0sin,,000,1,0,0,2uuxttxuxxuxtututt(1)1．在第三部分写出问题（1）三层显格式。2．根据你写出的差分格式，编写有限差分法程序。将所写程序放到第四部分。3．取0.1,0.1hh，分别将0.5,1.0,1.5,2.0t时刻的数值解画图显示。4.该问题的解析解为,cossinuxttx，将四个时刻的数值解的误差画图显示，对数值结果进行简单的讨论。三．实验原理、方法（算法）、步骤1、三层显格式建立由于题中0.1,0.1hh，0,1,0,2xt，取10,200NM，故令网比0.1rh，,0,1,2,10jxjhj，,0,1,,200ktkk，在,jkxt内网个点处，利用二阶中心差商得到如下格式：31111222222kkkkkkjjjjjjuuuuuuoohh（2）略去误差项得到：122211121kkkkkjjjjjurururuu（3）其中1,2,9,1,2,,199jk，局部截断误差为22oh。对于初始条件,0sinuxx，建立差分格式为：0sinsin,0,1,2,10jjuxjhj（4）对于初始条件,00uxt，利用中心差商，建立差分格式为：11110,=0,1,2,102jjjjuuuuj即，（5）对于边界条件0,1,0,0,2ututt，建立差分格式为：00,0,1,,200kkNuuk（6）将差分格式延拓使0k为内点，代入（3）得到的式子再与（5）联立消去1ju后整理得到：12020201111122jjjjurururu（7）综上（3）、（4）、（6）、（7）得到三层显格式如下：（局部截断误差为22oh）12221110120202011021,1,2,9,1,2,,199sinsin,0,1,2,10111,1,2,9220,0,1,,200kkkkkjjjjjjjjjjjkkNurururuujkuxjhjurururujuuk（8）其中0.1rh。四．实验环境（所用软件、硬件等）及实验数据文件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);end4%定义x=0边界上的数值fork=1:M+1u(k,1)=0;end%定义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,:)=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);%划分网格%作出数值解的函数图像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)精确解的函数图像');5xlabel('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*-');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;plot(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);%title('u(x,t)绝对误差的函数图像');%xlabel('x变量');%ylabel('t变量');%zlabel('绝对误差值');6五．实验结果及实例分析1、u(x,t)在t=0.5,1.0,1.5,2.0时刻的数值解、精确解以及绝对误差表1u(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.00200.00380.00520.00610.00640.00610.00520.00380.00200t=2.000.30900.58780.80900.95111.00000.95110.80900.58780.30900表2u(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.5878-0.30900t=1.500.00000.00000.00000.00000.00000.00000.00000.00000.00000t=2.000.30900.58780.80900.95111.00000.95110.80900.58780.30900表3u(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时刻的绝对误差均较小，由于0.11rh，该格式稳定，由数值计算得到的矩阵不难看出，数值解符合理论解。2、u(x,t)在t=0.5,1.0,1.5,2.0时刻的数值解、绝对误差函数图像图1数值解、精确解以及绝对误差函数图像7说明：上两图为函数的数值解与精确解，下两图为t=0.5,1.0,1.5,2.0时刻的数值解、绝对误差函数图像，符合理论解。教师签名年月日