《微分方程数值解》实验指导书2010

微分方程数值解实验指导书李光云数学与计算科学学院二O一O 十月一、概述本课程实验指导书是根据陆金甫，关治编著的《偏微分方程数值解法(第二版)》编写的。

通过上机实验，可帮助学生理解偏微分方程数值计算方法的基本思想，巩固学生所学过的算法，同时让学生进一步学习和掌握Matlab软件在数值计算方面中的应用。

二、实验环境本书选择的实验环境是计算机以及软件Matlab(版本6.5以上)一套。

三、实验课时安排共7次实验，16课时，试验五为综合设计性试验4课时，其它实验2个课时。

四、实验要求上机完成实验指导书中所规定的内容，自行按实验指导书要求完成程序设计和调试，并提交每次实验的实验报告，附带算法程序清单和算法输出结果。

五、实验考核要求上机完成试验内容，并提交一份算法程序清单和数值结果。

实验一 双曲型方程的迎风格式一、实验目的1、掌握双曲型方程的迎风格式的算法思想,格式稳定的条件；2、培养编程与上机调试能力。

二、实验课时：2个课时 三、实验准备1、了解偏微分方程的分类，以及双曲型方程的特点；2、熟悉求解偏微分方程有限差分法的基本原理，双曲型方程迎风格式的推导和特点，熟悉用迎风格式求解双曲型方程的算法步骤；3、熟悉Matlab 软件的基本命令及其数值计算中的基本命令和函数。

4、了解双曲型方程迎风格式稳定的条件。

四、实验内容用双曲型方程的迎风格式11110, 00, <0n n n nj jj j n n n njjj ju u u u a a h uuu u aa hττ+-++--+=>--+=其中，,h τ分别为时间步长和空间步长。

求解初值问题[]()()00,,0,,,0,u ux t T t xu x u x ∂∂⎧+=∈∈⎪∂∂⎨⎪=⎩其中()01,0,0,0.x u x x ≤⎧=⎨>⎩取10.01,2h λ==.计算[]0,1区间的u 值至0.5n t =,并画出解的图形。

五、基本思想及主要步骤考虑对流方程0,,0u u a x t t x∂∂+=∈>∂∂ 其a 为给定常数。

其迎风格式的基本思想是在方程中关于空间的偏导数用在特征线方向一侧的单边差商代替，得到对流方程的迎风格式11110, 00, <0n n n nj jj j n nn njjj ju u u u a a h uuu u aa hττ+-++--+=>--+=用Fourier 方法讨论格式的稳定性可知，以上两个格式都是条件稳定的，而且都在1,a hτλλ≤=时稳定。

程序设计步骤如下：(1)首先由a 的正负确定用哪个格式，再把格式改写为显示格式； (2)利用空间步长定义一个零向量用于存储u 值； (3)将初值离散存于u 向量内；(4)利用差分格式沿着t 坐标逐层计算u 值，每计算一次覆盖存储于u 向量内，直到算到0.5n t =层，然后输出结果；(5)利用得到的u 向量画出在0.5n t =时u 的函数图像。

六、实验提示：此算法可以编写M 程序，使得该程序可以实现对不同的a ，时间步长和空间步长求解，由于不用符号的a 用到的格式不同，迎风格式在条件1,a hτλλ≤=时稳定，故可以在程序里加一个判断部分，当不稳定时及时报警跳出。

七、实验结果记录和要求1、实验结果输出，画出解的图形；2、把得到的解和初值问题的解析解比较，观察比较结果；3、递交实验报告。

八、思考：若考虑对流方程的差分格式11110, <00, >0n n n nj jj j n n n njjj ju u u u a a h uuu u aa hττ+-++--+=--+=是否也能得到差分格式稳定？实验二 常系数扩散方程的经典差分格式一、实验目的1、了解抛物型方程的经典差分格式,格式稳定的条件；2、掌握常系数扩散方程初边值问题的加权隐式格式的求解方法；3、培养编程与上机调试能力。

二、实验课时：2个课时 三、实验准备1、了解偏微分方程的分类方法，抛物型方程的特点；2、了解常见的抛物型方程的差分格式；3、熟悉抛物型方程的加权隐式格式及其推导过程；4、熟悉常系数抛物型方程初边值问题的初值和边界离散方法。

四、实验内容考虑常系数扩散方程的初边值问题()()()22,01,0,,0sin ,01,0,1,0,0.u ux t t x u x x x u t u t t π⎧∂∂=<<>⎪∂∂⎪⎪=≤≤⎨⎪==≥⎪⎪⎩其解析解为()2,sin ,01,0.t u x t e x x t ππ-=≤≤≥用加权隐式格式近似求解()11111111222210n n n n n n n n j jj j j j j j u u u u u u u u h h θθτ----+-+-⎡⎤--+-+-++=⎢⎥⎢⎥⎣⎦其中01θ≤≤，取()110,,0,1,...,,10j J h x jh j J τ====为时间步长，2h τλ=为网格比，对不同的时间步长(11,,1,2,4,842λ=)，计算当10,1,2θ=时初边值问题的解()0.4,0.4u ，并且与精确解比较，分析比较结果。

五、基本思想及主要步骤用有限差分法解常系数扩散方程22u ua t x∂∂=∂∂有加权隐式差分格式()11111111222210n n n n n n n n j jj j j j j j u u u u u u u u a h h θθτ----+-+-⎡⎤--+-+-++=⎢⎥⎢⎥⎣⎦其中01θ≤≤，当12θ=时为Crank-Nicolson 格式，当1θ=时为向后差分格式，当0θ=时为向前差分格式。

加权隐式格式稳定的条件是12,121, 1.2a λθθθ≤≤-≤≤1，当0<2无限制当加权隐式格式是两层隐式格式，用第n 层计算第n+1层节点值的时候，要解线性方程组。

实验步骤如下：(1)输入.θλ，确定加权隐式格式的参数；(2)定义向量v ，把初边值条件离散，得到0,0,1,...,j u j J =的值存入向量v ； (3)利用差分格式由第n 层计算第n+1层，建立相应线性方程组，求解并且存入向量v ；(4)计算到0.4t =，输出()0.4,0.4u 。

六、实验提示：求解程序可编写M 文件，以便改变.θλ，得到不用参数下()0.4,0.4u 的值，然后比较与真实值的误差，判断格式的稳定性。

七、实验结果记录和要求1、输出不同参数下求得的节点的近似值；2、比较近似值和真实值的误差，印证格式稳定条件；3、实验结束后，递交实验报告。

八、思考：抛物型方程和双曲型方程的初边值问题的边界处理有什么不同，用隐式格式计算与显示格式相比，有什么优劣。

实验三 椭圆型方程的五点格式一、实验目的1、掌握poisson 方程定解问题的五点格式的求解方法；2、培养编程与上机调试能力。

二、实验课时：2个课时 三、实验准备1、了解偏微分方程的分类方法，椭圆型方程的特点；2、了解常见的poisson 方程的差分格式；3、熟悉poisson 方程的五点格式及其推导过程。

四、实验内容给定如下Laplace 方程（poisson 方程的特殊情况）的定解问题{}()()()()0,017,010,0,17,,00,,10100.u x y u y u y u x u x ∆=Ω=<<<<====利用椭圆型方程的五点格式，计算该问题的近似解，并且画出近似解的图形。

五、基本思想及主要步骤对Laplace 方程的第一边值问题()()()()22220,,,,,,u uu x y Dx y u x y x y x y Dα∂∂∆=+=∈∂∂=∈∂ 利用taylor 展开可得逼近它的五点差分格式的差分逼近()()1, 1.,1.122220,,,,i j ij i ji j ij i j h i j hij ij i j hu u u u u u u x y D h k u x y D α+-+--+-+=+=∈=∈∂其中,h k 分别为x 轴和y 轴步长，边界条件可以由()(),,u x y x y α=离散可得，当(,)x y D ∈∂时有(),ij i j x y αα=。

注意五点格式计算节点是由边界的已知节点，计算内部节点，计算时需要联立大型方程组，该方程组可以用迭代法求解。

主要步骤：(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 工具箱在求解其他偏微分方程定解问题的使用方法。

四、实验内容(1)一块有矩形裂纹的金属块,左侧被加热到100度，右侧以恒定速率降低到周围空气温度，所有其他边界独立，则有如下定解问题：0,,100,(Dirichlet ),10,(Neumann ),0,(Neumann )udu u tu un un∂-∆=∈Ω∂=∂=-∂∂=∂右侧条件左侧条件其它边界条件 设起始时刻0t 时刻金属块的温度为零度。

金属块长度为1.6，宽度为1，裂纹长度为0.8，宽度为0.1，指定起始时间为0，希望研究开始5s 内的热散发问题。

(2)一个内部充满空气的方框，内边长为0.2m ，外边长为0.5m 。

内边界上电势为1000v ，外边界上电势为0v ，域内没有电荷，现在求方框中的电势。

此问题即求方程0v ∇=边界条件为Dirichlet 条件，内边界上1000v =，外边界上0v =。

五、主要方法及提示这些问题都可归为求解偏微分方程的定解问题，我们用matlab 中的pde 工具箱求解此类问题的近似解。

对问题1，首先打开图形用户界面，画CSG 模型：首先画矩形（R1）四角坐标为[]0.5,0.5,0.5,0.5,x =--[]y 0.8,0.8,0.8,0.8=--，后画矩形（R2）四角坐标为[]=--，选择Options选项，打y0.4,0.4,0.4,0.4x=--[]0.05,0.05,0.05,0.05,开Grid Spacing对话框，在-0.05和0.05处输入x轴的附加短线，以帮助画出代表裂纹的矩形，然后显示网格和“Snap-to-gird”风格，金属块的CSG模型现在可以用公式R1-R2来表示。