抛物型方程数值解
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数学与计算科学学院 实 验 报 告
实验项目名称 所属课程名称 实 验 类 型 实 验 日 期 抛物型微分方程数值解 微分方程数值解 验证性 2015-4-23
班 学 姓 成
级 号 名 绩
信计 12-2 班 201253100212 黄全林
Biblioteka Baidu
一、实验概述： 【实验目的】
掌握抛物型方程的有限差分法，并学会应用有限差分法求解抛物型方程数值解 的 MATLAB 实现。本文利用向前差分法求解抛物型方程的数值解。
数值解与精确解在 x 方向上的走势趋同。
第 10 页 共 17 页
B. 空间 x 固定，观察时间 t 与 u 的走势，结果如下： 数值解：
精确解：
第 11 页 共 17 页
误差：
可以看出在时间方向基本趋同。
第 12 页 共 17 页
【实验小结】（收获体会） 通过此次上机亲自体验，使我对于数值计算有了更加深刻的认识。同时，我也感 觉自己的编程能力还有待提高，这对于数值计算是至关重要。这直接影响到了程序 的效率以及程序的运行时间以及计算精度。在以后的学习中，更应该加强这方面的 锻炼。 三、指导教师评语及成绩： 评语等级 评 语 优 良 中 及 格 不及格
n n u1n 2u0 u n n 1 sin t , u u , n 1 1 h2 n 1 n n n n uJ uJ uJ n n 1 2u J u J 1 sin tn , u J 1 u J 1 . 2 h
u x j ,0 cos x j , j 0,1, , J .
【实验环境】 1.硬件环境
第 1 页 共 17 页
2.2.软件环境
MATLAB7.0 二、实验内容： 【实验过程】（实验步骤） 实验任务 求解一维抛物方程的初边值问题：
u 2u sin t , 0 x 1, t 0, x x 2 u x 0, t u x 1, t 0, t 0, u x, 0 cos x, 0 x 1.
精确解：
第 4 页 共 17 页
误差：
数值解与精确解在 x 方向上的走势趋同。 B. 空间 x 固定，观察时间 t 与 u 的走势，结果如下：
第 5 页 共 17 页
数值解：
精确解：
误差：
第 6 页 共 17 页
可以看出在时间方向基本趋同。 （2）缩小步长，观察对输出的影响，取空间步长为 0.1，时间步长为 0.001，输 出结果为： 求解曲面：
第 16 页 共 17 页
hql_paowufangcheng.m 文件：
clc clear [p u e x t]=pwxywxq(0.1,0.005,10,200); figure
第 14 页 共 17 页
shading interp; surf(x,t,u); xlabel('x'); ylabel('t'); zlabel('u'); title('一维抛物线方程 向前欧拉法 数值解'); figure surf(x,t,p); xlabel('x'); ylabel('t'); zlabel('p'); title('精确解曲面'); figure surf(x,t,e); xlabel('x'); ylabel('t'); zlabel('e'); title('误差曲面'); figure plot(x,u) xlabel('x'); ylabel('u'); title('固定时间 figure plot(x,p) xlabel('x'); ylabel('p'); title('固定时间 figure plot(x,e) xlabel('x'); ylabel('u'); title('固定时间 figure plot(t,u) xlabel('t'); ylabel('u'); title('固定 x figure plot(t,p) 改变 t u与t 的关系 数值解图形') 改变 x u与x 的关系 误差图形') 改变 x u与x 的关系 精确解') 改变 x u与x 的关系 数值解')
精确解：u e
2t
cos x 1 cos t .
利用向前差分法求解 利用 MATLAB 进行求解，编辑函数文件 hql_xiangqianchafen.m。源程序见附录。 编辑调用函数 hql_xiangqianchafen( h,m,n,kmax,ep )的脚本文件，并作出相 应的求解曲面、精确解曲面和误差曲面图形。hql_paowufangcheng.m 源程序见附录。
第 15 页 共 17 页
xlabel('t'); ylabel('p'); title('固定 x figure plot(t,e) xlabel('t'); ylabel('e'); title('固定 x 改变 t u与t 的关系 误差图形')差曲面'); 改变 t u与t 的关系 精确解图形')
第 7 页 共 17 页
精确解曲面：
误差曲面：
第 8 页 共 17 页
通过相应曲面的对比，可以发现，此时精确解与数值解图形逼近效果明显趋同。 下面分别固定空间 x 和时间 t，观察数值解与精确解的逼近程度。 B. 时间固定，空间 x 与 u 的关系： 数值解：
第 9 页 共 17 页
精确解：
误差：
1.实验报告按时完成,字迹清楚,文字叙述流畅,逻辑性强 2.实验方案设计合理 3.实验过程（实验步骤详细,记录完整,数据合理,分析透彻） 4 实验结论正确.
成
绩： 指导教师签名： 批阅日期：
附录：源 程 序
第 13 页 共 17 页
函数文件： 文件名：hql_xiangqianchafen.m
function [p u e x t]=hql_xiangqianchafen(h1,h2,m,n) （Ut-aUxx=f(x,t),a>0) % 解抛物线型一维方程 向前欧拉格式
【实验原理】 抛物型方程的先前差分法求解： 方程离散：
1 un un j j
边值条件：
n n un j 1 2u j u j 1
h2
sin t n , j 1, , J 1, n 1, , N
j 0，
n 1 n u0 u0
j J,
初值条件：
%不用解线性方程组，由下一层(时间层）的值就直接得到上一层的值 %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); 定 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 %此处 r=a*h2/(h1*h1）；a=1 要求 r<=1/2 差分格式才稳
运行上述 hql_paowufangcheng.m 文件： （1）当空间步长为 0.1，时间步长为 0.005 时，结果如下：
第 2 页 共 17 页
求解曲面：
精确解曲面：
误差曲面：
第 3 页 共 17 页
通过相应曲面的对比，可以发现，此时精确解与数值解图形逼近效果明显趋同。 下面分别固定空间 x 和时间 t，观察数值解与精确解的逼近程度。 A. 时间固定，空间 x 与 u 的关系： 数值解：
实验项目名称 所属课程名称 实 验 类 型 实 验 日 期 抛物型微分方程数值解 微分方程数值解 验证性 2015-4-23
班 学 姓 成
级 号 名 绩
信计 12-2 班 201253100212 黄全林
Biblioteka Baidu
一、实验概述： 【实验目的】
掌握抛物型方程的有限差分法，并学会应用有限差分法求解抛物型方程数值解 的 MATLAB 实现。本文利用向前差分法求解抛物型方程的数值解。
数值解与精确解在 x 方向上的走势趋同。
第 10 页 共 17 页
B. 空间 x 固定，观察时间 t 与 u 的走势，结果如下： 数值解：
精确解：
第 11 页 共 17 页
误差：
可以看出在时间方向基本趋同。
第 12 页 共 17 页
【实验小结】（收获体会） 通过此次上机亲自体验，使我对于数值计算有了更加深刻的认识。同时，我也感 觉自己的编程能力还有待提高，这对于数值计算是至关重要。这直接影响到了程序 的效率以及程序的运行时间以及计算精度。在以后的学习中，更应该加强这方面的 锻炼。 三、指导教师评语及成绩： 评语等级 评 语 优 良 中 及 格 不及格
n n u1n 2u0 u n n 1 sin t , u u , n 1 1 h2 n 1 n n n n uJ uJ uJ n n 1 2u J u J 1 sin tn , u J 1 u J 1 . 2 h
u x j ,0 cos x j , j 0,1, , J .
【实验环境】 1.硬件环境
第 1 页 共 17 页
2.2.软件环境
MATLAB7.0 二、实验内容： 【实验过程】（实验步骤） 实验任务 求解一维抛物方程的初边值问题：
u 2u sin t , 0 x 1, t 0, x x 2 u x 0, t u x 1, t 0, t 0, u x, 0 cos x, 0 x 1.
精确解：
第 4 页 共 17 页
误差：
数值解与精确解在 x 方向上的走势趋同。 B. 空间 x 固定，观察时间 t 与 u 的走势，结果如下：
第 5 页 共 17 页
数值解：
精确解：
误差：
第 6 页 共 17 页
可以看出在时间方向基本趋同。 （2）缩小步长，观察对输出的影响，取空间步长为 0.1，时间步长为 0.001，输 出结果为： 求解曲面：
第 16 页 共 17 页
hql_paowufangcheng.m 文件：
clc clear [p u e x t]=pwxywxq(0.1,0.005,10,200); figure
第 14 页 共 17 页
shading interp; surf(x,t,u); xlabel('x'); ylabel('t'); zlabel('u'); title('一维抛物线方程 向前欧拉法 数值解'); figure surf(x,t,p); xlabel('x'); ylabel('t'); zlabel('p'); title('精确解曲面'); figure surf(x,t,e); xlabel('x'); ylabel('t'); zlabel('e'); title('误差曲面'); figure plot(x,u) xlabel('x'); ylabel('u'); title('固定时间 figure plot(x,p) xlabel('x'); ylabel('p'); title('固定时间 figure plot(x,e) xlabel('x'); ylabel('u'); title('固定时间 figure plot(t,u) xlabel('t'); ylabel('u'); title('固定 x figure plot(t,p) 改变 t u与t 的关系 数值解图形') 改变 x u与x 的关系 误差图形') 改变 x u与x 的关系 精确解') 改变 x u与x 的关系 数值解')
精确解：u e
2t
cos x 1 cos t .
利用向前差分法求解 利用 MATLAB 进行求解，编辑函数文件 hql_xiangqianchafen.m。源程序见附录。 编辑调用函数 hql_xiangqianchafen( h,m,n,kmax,ep )的脚本文件，并作出相 应的求解曲面、精确解曲面和误差曲面图形。hql_paowufangcheng.m 源程序见附录。
第 15 页 共 17 页
xlabel('t'); ylabel('p'); title('固定 x figure plot(t,e) xlabel('t'); ylabel('e'); title('固定 x 改变 t u与t 的关系 误差图形')差曲面'); 改变 t u与t 的关系 精确解图形')
第 7 页 共 17 页
精确解曲面：
误差曲面：
第 8 页 共 17 页
通过相应曲面的对比，可以发现，此时精确解与数值解图形逼近效果明显趋同。 下面分别固定空间 x 和时间 t，观察数值解与精确解的逼近程度。 B. 时间固定，空间 x 与 u 的关系： 数值解：
第 9 页 共 17 页
精确解：
误差：
1.实验报告按时完成,字迹清楚,文字叙述流畅,逻辑性强 2.实验方案设计合理 3.实验过程（实验步骤详细,记录完整,数据合理,分析透彻） 4 实验结论正确.
成
绩： 指导教师签名： 批阅日期：
附录：源 程 序
第 13 页 共 17 页
函数文件： 文件名：hql_xiangqianchafen.m
function [p u e x t]=hql_xiangqianchafen(h1,h2,m,n) （Ut-aUxx=f(x,t),a>0) % 解抛物线型一维方程 向前欧拉格式
【实验原理】 抛物型方程的先前差分法求解： 方程离散：
1 un un j j
边值条件：
n n un j 1 2u j u j 1
h2
sin t n , j 1, , J 1, n 1, , N
j 0，
n 1 n u0 u0
j J,
初值条件：
%不用解线性方程组，由下一层(时间层）的值就直接得到上一层的值 %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); 定 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 %此处 r=a*h2/(h1*h1）；a=1 要求 r<=1/2 差分格式才稳
运行上述 hql_paowufangcheng.m 文件： （1）当空间步长为 0.1，时间步长为 0.005 时，结果如下：
第 2 页 共 17 页
求解曲面：
精确解曲面：
误差曲面：
第 3 页 共 17 页
通过相应曲面的对比，可以发现，此时精确解与数值解图形逼近效果明显趋同。 下面分别固定空间 x 和时间 t，观察数值解与精确解的逼近程度。 A. 时间固定，空间 x 与 u 的关系： 数值解：