《微分方程数值解》实验指导书.docx
《《微分方程数值解》实验指导书.docx》由会员分享，可在线阅读，更多相关《《微分方程数值解》实验指导书.docx（22页珍藏版）》请在冰点文库上搜索。
《微分方程数值解》实验指导书
微分方程数值解
实验指导书
李光云
数学与计算科学学院
O一0十月
一、概述
本课程实验指导书是根据陆金甫，关治编著的《偏微分方程数值解法（第二版）》编写的。
通过上机实验，可帮助学生理解偏微分方程数值计算方法的基本思想，巩固学生所学过的算法，同时让学生进一步学习和掌握Matlab软件在数值计算方面中的应用。
二、实验环境
本书选择的实验环境是计算机以及软件Matlab（版本6.5以上）一套。
三、实验课时安排
共7次实验，16课时，试验五为综合设计性试验4课时，其它实验2个课时。
四、实验要求
上机完成实验指导书中所规定的内容，自行按实验指导书要求完成程序设计和调试，并提交每次实验的实验报告，附带算法程序清单和算法输出结果。
五、实验考核要求
上机完成试验内容，并提交一份算法程序清单和数值结果。
实验一双曲型方程的迎风格式
、实验目的
1、掌握双曲型方程的迎风格式的算法思想，格式稳定的条件;
2、培养编程与上机调试能力。
、实验课时：
2个课时
三、实验准备
1、了解偏微分方程的分类，以及双曲型方程的特点；
2、熟悉求解偏微分方程有限差分法的基本原理，双曲型方程迎风格式的推导和特点，熟悉用迎风格式求解双曲型方程的算法步骤；
3、熟悉Matlab软件的基本命令及其数值计算中的基本命令和函数。
4、了解双曲型方程迎风格式稳定的条件。
四、实验内容
用双曲型方程的迎风格式
n1n
UjUj
nn
a0
a<0
Uj-Uj」
0,h
n1nnn
Uj—UjUj1—Ujc
a0,
th
其中，•,h分别为时间步长和空间步长。
求解初值问题
賦”,ts,
Iu（x,0）=u°（x）,
其中
p!
x",
U。
八0,x0.
五、基本思想及主要步骤
考虑对流方程
:
U:
u
a0,x.,t0
**
-t:
x
其a为给定常数。
其迎风格式的基本思想是在方程中关于空间的偏导数用在特征线
方向一侧的单边差商代替，得到对流方程的迎风格式
n1nnn
Uj_UjaUj_U2二0,
n1nnn
Uj~UjaUj1_出=0,
a0
a<0
用FoUrier方法讨论格式的稳定性可知，以上两个格式都是条件稳定的，而且
都在a九兰1,九=—时稳定。
h
程序设计步骤如下：
（1）首先由a的正负确定用哪个格式，再把格式改写为显示格式；
⑵利用空间步长定义一个零向量用于存储u值；
（3）将初值离散存于U向量内；
（4）利用差分格式沿着t坐标逐层计算U值，每计算一次覆盖存储于U向量内,
直到算到tn=0.5层，然后输出结果;
（5）利用得到的U向量画出在tn=0.5时U的函数图像。
六、实验提示：
此算法可以编写M程序，使得该程序可以实现对不同的a，时间步长和空间
步长求解，由于不用符号的a用到的格式不同，迎风格式在条件a九兰1花=上时
h
稳定，故可以在程序里加一个判断部分，当不稳定时及时报警跳出。
七、实验结果记录和要求
1、实验结果输出，画出解的图形；
2、把得到的解和初值问题的解析解比较，观察比较结果;
3、递交实验报告。
八、思考：
若考虑对流方程的差分格式
n“1n
Uj-Uj
T
n1n
uUj
T
nn
Uj—Uj—
0,h
nn
Uj1—Ujc
a0,
h
a<0
a>0
是否也能得到差分格式稳定?
实验二常系数扩散方程的经典差分格式
、实验目的
1、了解抛物型方程的经典差分格式，格式稳定的条件；
2、掌握常系数扩散方程初边值问题的加权隐式格式的求解方法;
3、培养编程与上机调试能力。
、实验课时：
2个课时
三、实验准备
1、了解偏微分方程的分类方法，抛物型方程的特点；
2、了解常见的抛物型方程的差分格式；
3、熟悉抛物型方程的加权隐式格式及其推导过程；
4、熟悉常系数抛物型方程初边值问题的初值和边界离散方法。
四、实验内容
考虑常系数扩散方程的初边值问题
L2
cucu
——=一,0VX<1,tA0,猊cX
u（x,0）=sin兀x,0兰x兰1,
u0,t二u1,t=0,t_0.
其解析解为
用加权隐式格式近似求解
1T
其中0一二一1,取J=10,h,Xj二jhj=0,1,...,J,•为时间步长，2为
10h
111
网格比，对不同的时间步长（，丄，1,2,4,8），计算当^-0,1」时初边值问题的
422
解u0.4,0.4，并且与精确解比较，分析比较结果。
五、基本思想及主要步骤
■■-2
U
用有限差分法解常系数扩散方程
：
:
u
；:
t
有加权隐式差分格式
1
其中0冬二乞1,当时为Crank-Nicolson格式，当v-1时为向后差分格式，当
2
v-0时为向前差分格式。
加权隐式格式稳定的条件是
11
2a，当0乞T,
1-2日2
无限制，当1"<1.
2
加权隐式格式是两层隐式格式，用第n层计算第n+1层节点值的时候，要解线
性方程组。
实验步骤如下：
（1）输入='，确定加权隐式格式的参数；
（2）定义向量v，把初边值条件离散，得到U0,j=0,1,...,J的值存入向量V；
（3）利用差分格式由第n层计算第n+1层，建立相应线性方程组，求解并且存入
向量v；
（4）计算到t=0.4，输出u0.4,0.4。
六、实验提示：
求解程序可编写M文件，以便改变=■,得到不用参数下u0.4,0.4的值,然后比较与真实值的误差，判断格式的稳定性。
七、实验结果记录和要求
1、输出不同参数下求得的节点的近似值；
2、比较近似值和真实值的误差，印证格式稳定条件；
3、实验结束后，递交实验报告。
八、思考：
抛物型方程和双曲型方程的初边值问题的边界处理有什么不同，用隐式格式计算与显示格式相比，有什么优劣。
实验三椭圆型方程的五点格式
、实验目的
1、掌握poisson方程定解问题的五点格式的求解方法;
2、培养编程与上机调试能力。
、实验课时：
2个课时
三、实验准备
1、了解偏微分方程的分类方法，椭圆型方程的特点;
2、了解常见的poisson方程的差分格式；
3、熟悉poisson方程的五点格式及其推导过程。
四、实验内容
给定如下
Laplace方程（poisson方程的特殊情况）的定解问题u=0,门-10：
：
x:
:
17,0：
：
y：
10二
u0,y]=u17,yi=ux,0j=0,
ux,10=100.
利用椭圆型方程的五点格式，计算该问题的近似解，并且画出近似解的图形。
五、基本思想及主要步骤
对Laplace方程的第一边值问题
ux,y二[x,y,x,y
利用taylor展开可得逼近它的五点差分格式的差分逼近
其中h,k分别为x轴和y轴步长，边界条件可以由ux,yx,y离散可得，当
（x,y）VD时有九="[Xi,yj。
注意五点格式计算节点是由边界的已知节点，计算内部节点，计算时需要联立大型方程组，该方程组可以用迭代法求解。
主要步骤：
（1）首先取定h,k，对求解区域划分网格，按照网格定义矩阵v1，使得矩阵里面
每个元素对应求解区域中的每个节点；
（2）由边界条件定义v1矩阵中边界元素的值，其余元素定义为零；
（3）定义与v1同型的零矩阵v2；
（4）用五点格式公式通过矩阵v1迭代计算矩阵v2,迭代精度为0.1；
（5）画图。
六、实验提示：
由v1迭代计算矩阵v2时，v2元素通过v1元素利用五点格式计算，当v1v2对
应元素误差在控制范围之内时，终止迭代，否则继续，直到满足误差条件。
七、实验结果记录和要求
1、输出所有节点的计算结果，并且画图；
2、实验结束后，递交实验报告。
八、思考：
如果求解区域是一般区域，如何划分网格，处理边界条件。
实验四pde工具箱在热传导，静电学等方面的应用
一、实验目的
1、全面了解pde工具箱对其他偏微分方程的求解方法；
2、培养编程与上机调试能力。
二、实验课时：
2个课时
三、实验准备
1、了解静电学，热传导等方面的偏微分方程；
2、熟悉pde工具箱的参数的具体意义；
3、熟悉Matlab软件中pde工具箱在求解其他偏微分方程定解问题的使用方法。
四、实验内容
⑴一块有矩形裂纹的金属块，左侧被加热到100度，右侧以恒定速率降低到周围空气温度，所有其他边界独立，则有如下定解问题：
:
u
du=0,u■」，
；:
t
u=100,右侧（Dirichlet条件）,
:
u
'10,左侧（Neumanr条件）,
fn
:
u
-=0,其它边界（Neumann条件）
in
1.6，宽度为1,裂纹长度为
设起始时刻t0时刻金属块的温度为零度。
金属块长度为
0.8，宽度为0.1，指定起始时间为0,希望研究开始5s内的热散发问题。
⑵一个内部充满空气的方框，内边长为0.2m，外边长为0.5m。
内边界上电势
为1000v，外边界上电势为0v，域内没有电荷，现在求方框中的电势。
此问题即求方程
'v=0
边界条件为Dirichlet条件，内边界上v=1000，外边界上v=0。
五、主要方法及提示
这些问题都可归为求解偏微分方程的定解问题，我们用matlab中的pde工具箱
求解此类问题的近似解。
对问题1,首先打开图形用户界面，画CSG模型：
首先画矩形（R1）四角坐标为x-丨-0.5,0.5,0.5,-0.51,
y-丨-0.8,-0.8,0.8,0.8丨，后画矩形（R2）四角坐标为
x-丨-0.05,0.05,0.05,-0.051,y-丨-0.4,-0.4,0.4,0.4丨，选择Options选项，打
开GridSpacing对话框，在-0.05和0.05处输入x轴的附加短线，以帮助画出代表裂纹的矩形，然后显示网格和“Snap-to-gird”风格，金属块的CSG模型现在可以用公式R1-R2来表示。
然后指定边界条件，输入pde参数，初始网格，SolveParameters对话框中输入
时间参数，在不同时间段求解方程。
可以显示动态传导过程，首先在PlotSelection对话框中选择Animation，选择colormaphot，单击Plot，在单独的图形窗口中开始解的图形记录。
问题2,首先打开用户图形界面，选择Electrostatics模式并建立几何模型，首先画边长为0.2的方形SQ1，然后画一个中心位置相同的边长为0.5的方形SQ2,
则分析域为SQ2-SQ1。
在Setformula框中输入表达式，输入边界条件，打开PDESpecification对话框，在Spacechargedensity中输入0.Dielectricity设为1.
初始化网格，求解方程。
六、实验结果记录和要求
1、记录检验结果，保存M文件，保存输出图像；
2、实验结束后，递交实验报告。
七、思考：
Matlab软件的Pde工具箱在求解偏微分方程中的作用与优势。
实验五线性椭圆型方程的有限元方法实验
、实验目的
1、掌握线性椭圆型方程有限元方法的基本思想；
2、掌握用matlab的pde工具箱求解线性椭圆型方程的方法；
3、掌握用命令行求解线性椭圆型方程的方法及其相关命令；
4、掌握用三角线性元求解线性椭圆型方程的MATLAB解法。
、实验课时：
4个课时
三、实验准备
1、熟悉线性椭圆型方程的一般形式；
2、掌握有限元方法的基本思想和推导过程以及三角形线性元求解的思路和方法；
3、熟悉Matlab软件的基本命令函数的用法；
4、掌握pde工具箱参数的设定和各个菜单的用法。
四、实验内容
1.用pde工具箱和命令行先后求解椭圆型方程边值问题
-ux,yi=0,0cu了兀x'
ux,0;=u0,y=5,yi=0,ux,10;=100sin
ex（10丿
该问题的精确解为
100shiysinix
ux,y=
l10丿J0丿
sh■:
比较求得的近似解和精确解的误差，并且改变网格粗细观察误差的变化情况。
2.用matlab语言求解上述边值问题。
进行三角形网格划分：
五、实验方法与提示：
1.用pde工具箱求解，在命令窗口输入pdetool，开启工具箱，选用Generic
Scalar模式，
按顺序先画出求解区域为单位圆，
设置边界条件，在BoundaryCondition对话框中选择Dirichler条件。
设置方程参数，在pdeSpecilication对话框中选择Elliptic方程。
戈U分网格，还可以加密网格。
求解，通过菜单画图，与精确解比较误差，输出网格节点编号，单元编号，节点坐标，以及解的数值。
初始网格求解完成后，还可以细分网格，在比较误差，看网格细分后，误差的变化。
2.用命令求解，提示所用：
函数circleg.m返回单位圆盘边界点坐标，circleb.m描述边界条件，命令initmesh（）初始划分网格，refinemesh（）细分网格，assempde（）求解方程，pdemesh（）生成网格图，pdesurf（）生成解表面图及误差表面图。
可以在Help菜单中得到函数的详细用法。
3.用matlab语言求解，需按照有限元的计算方法，首先划分网格，然后在每
个单元上计算单元刚度矩阵和单元荷载向量，再总体合成总刚度矩阵和总荷载向量，
再考虑边界条件，最后求解线性方程组。
六、实验结果记录和要求
1、生成程序清单，输出近似解和误差的图像；
2、用pde工具箱求解保存M文件；
3、实验结束后，递交实验报告。
七、思考：
有限元方法求解偏微分方程和有限差分法求解有什么不同。
附件一
算法举例
问题一
考虑初值问题
Q讥0
tiex
ux,0i=Uox
其中
4xe[0.4,0.6]
UoX=0,x[0.4,0.6]
用迎风格式计算上述冋题。
（考虑空间步长h=0.01，网格比0.5，对不同
h
的t值计算u，并且画出u-X图像）
1
解:
该初值问题的迎风格式可写为简单形式：
U；1U；•U；」.
注意到我们只在非零区间计算U值。
用matlab编写程序：
function[u,t]=yf（k）
u=ones（1,22）;u
（1）=0;%0.4〜0.6值为1，空间步长为0.01
v=zeros（1,22）;
forn=1:
k
u（n+22）=0;v（n+22）=0;%定义1后为零的元素
fori=2:
n+22
v（i）=（u（i）+u（i-1））/2;%用迎风格式求解
i=i+1;
end
w=u;
u=v;
v=w;
n=n+1;
end
x=0.39:
0.01:
（0.6+k*0.01）;
t=k*0.005;
plot（x,u）;
在matlab命令:
>>[u,t]=yf（10）
输出结果为：
t=
0.0500
图像输出为:
在matlab命令：
>>[u,t]=yf（100）输出结果为：
t=
0.5000
图像输出为:
并且可以改变不同的k，得到不同时刻t的u值，进而可以观察u值随着时间t的变化而变化的趋势。
问题二
求单位圆盘泊松方程的边值问题
22
rux,y=1,xy：
：
1,
22
ux,y=0,xy1.
该问题的精确解为
1.用图形用户界面计算
在命令窗口中输入pdetool命令，选用GenericScalar模式
（1）单击Option菜单，选择addagrid选项，设置"snap-to-grid”特点。
单击快捷键画圆形区域，若该图形不是单位圆，双击该圆，打开对话框，在其中可以指定圆心的精确位置和半径。
（2）通过快捷键设置边界模式。
分割的几何边界显示出来，并且外边界指定
为默认设置，即Dirichlet边界条件，u=0。
本例采用默认设置，若边界条件不同，可以通过双击边界打开一对话框，在其中输入对应的边界条
件。
（3）单击快捷键，定义偏微分方程，该操作打开一对话框，可以在其中定义
PDE系数c,a和f。
本例中，它们均为常数：
c=1，f=1，a=0。
（4）选择Mesh菜单中的InitializeMesh选项，初始化显示三角形网格。
（5）在Mesh菜单中选择RefineMesh选项，改进初始网格并显示新网格。
（6）单击求解快捷键进行方程求解，Matlab可以用图形来表示问题的解，单
击图形快捷键，打开PlotSelection对话框。
利用该对话框，可以选择不同类型的解图。
⑺比较数值解和精确解之间的误差。
在PlotSelection对话框中的Color选
择Property弹出式菜单中的userentry。
然后在Userentry编辑区中输入精确解表达式。
可以获得解的绝对误差的图形表示。
Color：
uHeightu
以上为解的表面图形和误差图形
2.使用命令行函数求解在命令窗口输入：
>>[p,e,t]=initmesh（Circleg','Hmax',1）;
>>error=[];err=1;
>>whileerr>0.001,
[p,e,t]=refinemesh（Circleg',p,e,t）;
u=assempde（cirCleb1’,p,e,t,1,0,1）;
exact=-（p（1,:
）.A2+p（2,:
）.A2-1）/4;
err=norm（u-exact',inf）;
error=[error,err];
end
>>pdemesh（p,e,t）;%生成网格图
>>pdesurf（p,t,u）;%生成解的表面图
>>pdesurf（p,t,u-exact'）;%生成解的绝对误差的表面图
上面的命令行中，第一行用参数化函数circleg创建初始网格。
为解的最大误差设定初始向量error，设置初始误差err为1。
此循环将一直进
行下去，直到解的误差小于0.001。
（1）改进网格；
（2）求解线性系统，注意本题中椭圆型PDE的系数为常数（c=f=1,a=0）,circleb1包含边界条件的描述，p,e,t定义三角形网格；
（3）求数值解和精确解的误差。
Exact向量包含节点处的精确解；
（4）绘制网格，解和误差的图形。
附件二
实验报告的撰写要求
一、写出实验课程名称、课号，任课老师;
二、写出姓名、学号、日期；
三、写出实验目的、实验内容；
四、写出实验结果，分析实验结果；
五、写出程序；
六、写出心得体会。
升级会员 