抛物型方程数值解

数学与计算科学学院实验报告

实验项目名称抛物型微分方程数值解

所属课程名称微分方程数值解

实验类型验证性

实验日期2015-4-23

班级信计12-2班

学号201253100212

姓名黄全林

成绩

一、实验概述：

【实验目的】

掌握抛物型方程的有限差分法，并学会应用有限差分法求解抛物型方程数值解的MATLAB 实现。本文利用向前差分法求解抛物型方程的数值解。

【实验原理】

抛物型方程的先前差分法求解：

方程离散：

N

n J j t h u u u u u n n j n j n j n j

n j ,,1,1,,1,sin 22111⋅⋅⋅=-⋅⋅⋅=++-=--++τ

边值条件：0=j ，100

1011122sin ,,n n n n n n n n u u u u u t u u h τ+----+=+=,

J j =1111122sin ,.n n n n n n n J J J J J n J J u u u u u t u u h

τ++-+---+=+=初值条件：().

,,1,0,cos 0,J j x x u j j ⋅⋅⋅==π【实验环境】

1.硬件环境

2.2.软件环境

MATLAB7.0

二、实验内容：

【实验过程】（实验步骤）

实验任务

求解一维抛物方程的初边值问题：

()()()()()

222sin ,01,t 0,0,t 1,0,t 0,

x,0cos ,0 1.

cos 1cos .x x t u u t x x x

u u t u x x u e x t πππ-∂∂=+<<>∂∂==>=<<=+-精确解：利用向前差分法求解

利用MATLAB 进行求解，编辑函数文件hql_xiangqianchafen.m。源程序见附录。编辑调用函数hql_xiangqianchafen(h,m,n,kmax,ep )的脚本文件，并作出相应的求解曲面、精确解曲面和误差曲面图形。hql_paowufangcheng.m 源程序见附录。运行上述hql_paowufangcheng.m 文件：

（1）当空间步长为0.1，时间步长为0.005时，结果如下：

求解曲面：

精确解曲面：

误差曲面：

通过相应曲面的对比，可以发现，此时精确解与数值解图形逼近效果明显趋同。下面分别固定空间x和时间t，观察数值解与精确解的逼近程度。

A.时间固定，空间x与u的关系：

数值解：

精确解：

误差：

数值解与精确解在x方向上的走势趋同。

B.空间x固定，观察时间t与u的走势，结果如下：

数值解：

精确解：

误差：

可以看出在时间方向基本趋同。

（2）缩小步长，观察对输出的影响，取空间步长为0.1，时间步长为0.001，输出结果为：

求解曲面：

精确解曲面：

误差曲面：

通过相应曲面的对比，可以发现，此时精确解与数值解图形逼近效果明显趋同。下面分别固定空间x和时间t，观察数值解与精确解的逼近程度。

B.时间固定，空间x与u的关系：

数值解：

精确解：

误差：

数值解与精确解在x方向上的走势趋同。

B.空间x固定，观察时间t与u的走势，结果如下：数值解：

精确解：

误差：

可以看出在时间方向基本趋同。

【实验小结】（收获体会）

通过此次上机亲自体验，使我对于数值计算有了更加深刻的认识。同时，我也感觉自己的编程能力还有待提高，这对于数值计算是至关重要。这直接影响到了程序的效率以及程序的运行时间以及计算精度。在以后的学习中，更应该加强这方面的锻炼。

三、指导教师评语及成绩：

评语

评语等级

优良中

及

格

不及格

1.实验报告按时完成,字迹清楚,文字叙述流畅,逻辑性强

2.实验方案设计合理

3.实验过程（实验步骤详细,记录完整,数据合理,分析透彻）

4实验结论正确.

成绩：

指导教师签名：

批阅日期：附录：源程序

函数文件：

文件名：hql_xiangqianchafen.m

function[p u e x t]=hql_xiangqianchafen(h1,h2,m,n)

%解抛物线型一维方程向前欧拉格式（Ut-aUxx=f(x,t),a>0)

%不用解线性方程组，由下一层(时间层）的值就直接得到上一层的值

%m,n为x,t方向的网格数，例如（2-0）/0.01=200;

%e为误差，p为精确解

u=zeros(n+1,m+1);

x=0+(0:m)*h1;

t=0+(0:n)*h2;

for(i=1:n+1)

u(i,1)=0;

u(i,m+1)=0;

end

for(i=1:m+1)

u(1,i)=cos(pi*x(i));

end

for(i=1:n+1)

for(j=1:m+1)

f(i,j)=0;

end

end

r=h2/(h1*h1);%此处r=a*h2/(h1*h1）；a=1要求r<=1/2差分格式才稳定

for(i=1:n)

for(j=2:m)

u(i+1,j)=(1-2*r)*u(i,j)+r*(u(i,j-1)+u(i,j+1))+h2*f(i,j);

end

end

for(i=1:n+1)

for(j=1:m+1)

p(i,j)=exp(-pi*pi*t(i))*cos(pi*x(j))+(1-cos(t(i)));

e(i,j)=abs(u(i,j)-p(i,j));

end

end

hql_paowufangcheng.m文件：

clc

clear

[p u e x t]=pwxywxq(0.1,0.005,10,200);

figure