实验四 微分方程符号解与数值解

实验四、 微分方程的符号解与数值解

【实验目的】

1．求微分方程的符号解。

2．求微分方程的数值解。

【实验内容】

问题1：求4290y y y '''++=满足初始条件(0)0,(0)15y y '==的符号解。

问题2：求解微分方程22,(0)2,0104(1)

y t y y t t --'==≤≤+.先求符号解,再求数值解, 并作图进行比较.

问题3：求解微分方程

2(1)2,(0)1,(0)3x y xy y y ''''+===.先求符号解,再求数值解, 并进行比较.

【相关介绍】

(一)微分方程的符号解.

微分方程的符号解也叫做解析解. 求微分方程（组）的符号解用命令dsolve. 命令格式如下：

说明：用字符串表示方程，自变量缺省则默认为t. 导数用D 表示, 2阶导数用D2表示,以此类推。返回值s 是符号解.

例1: 求y ay b '=+的符号解.

s=dsolve('Dy=a*y+b') %自变量缺省则默认为t

例2: 求sin(2)y x y ''=-满足初始条件(0)0,(0)1y y '==的符号解.

s=dsolve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x')

simplify(s) %如果得到符号解比较复杂,可以试试化简 大家得到的结果是什么呢？

(二)微分方程的数值解.

一般说来,只有对一些典型的常微分方程,才能求出它们的一般解.然而在实际问题中遇到的常微分方程往往很复杂,在许多情况下得不出一般解.所以一般是要求在若干点的近似数值解. 求数值解的命令如下:

说明：(1) 返回值中，xout 表示自变量的取值点(x 0,x 1,…,x n )', yout 表示数值解,它是一个矩阵,它的每一列对应y 的一个分量。

(2) 这里'equation'必须是事先定义的表示微分方程(组)的M-文件。

(3) [x0,xm]是自变量的区间。

(4) y0是初始向量值。

(5)ode45还可以换成其他算法,如ode23.

注意: 命令ode45或ode23是对一阶常微分方程(组)设计的,因此对高阶常微分方程,需将它转化为一阶常微分方程组.

例如对二阶常微分方程(1)0y y y y '''+++=,通过令

12,y y y y '==得到一阶常微分方方程组11221

(1)y y y y y y ?'=-+-??'=??.

在求数值解时,我们往往将数值解与画图结合, 将数值解用图像呈现出来.

例3: 教材《数学模型》第140页的数值计算. 实际上是求解微分方程组1

121

212

120.3(0)0.02,(0)0.98

dx x x x dt dx x x dt x x ?=-??

?=-??==???

先定义M-文件ill.m 。 function y=ill(t,x)

y=[x(1)*x(2)-0.3*x(1),-x(1)*x(2)]’;

然后

ts=[0,50]; % 课本是ts=0:50,也行

x0=[0.02,0.98];

[t,x]=ode45('ill',ts,x0);

plot(t,x(:,1),t,x(:,2)),

%根据x 的第一、二列同时作两条曲线

grid %为了观察方便,可添上网格线

例4: 求解微分方

1,(0)1y y t y '=-++=.先求符号解,再求数值解, 并作图

进行比较.

s=dsolve('Dy=-y+t+1','y(0)=1','t')

simplify(s)

可得符号解为y=t+exp(-t).

为了求数值解,先编写M-文件fun.m

function f=fun(t,y)

f=-y+t+1;

保存，再运行如下命令:

clear;close;t=0:0.02:1;

y=t+exp(-t);

plot(t,y) %画符号解的图形

hold on %保留已画好的图形,

[t,y]=ode45('fun',[0,1],1);

plot(t,y,'ro'); %画数值解图形,用红色小圆圈

xlabel('t'),ylabel('y') %标上各坐标名称

运行结果见下图, 可见符号解和数值解吻合得很好.

微分方程数值解实验报告

微分方程数值解法 课程设计报告 班级：_______ 姓名：___ 学号：__________ 成绩： 2017年 6月 21 日

摘要 自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge—Kutta方法、Adams法以及椭圆型方程、抛物型方程的有限差分方法等，通过具体的算例，结合MATLAB求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。 关键词：微分方程数值解、MATLAB 目录

摘要 (2) 目录 (3) 第一章常微分方程数值解法的基本思想与原理 (4) 1.1常微分方程数值解法的基本思路 (4) 1.2用matlab编写源程序 (4) 1.3常微分方程数值解法应用举例及结果 (5) 第二章常系数扩散方程的经典差分格式的基本思想与原理 (6) 2.1常系数扩散方程的经典差分格式的基本思路 (6) 2.2 用matlab编写源程序 (7) 2.3常系数扩散方程的经典差分格式的应用举例及结果 (8) 第三章椭圆型方程的五点差分格式的基本思想与原理 (10) 3.1椭圆型方程的五点差分格式的基本思路 (10) 3.2 用matlab编写源程序 (10) 3.3椭圆型方程的五点差分格式的应用举例及结果 (12) 第四章总结 (12) 参考文献 (12)

偏微分方程数值解实验报告

偏微分方程数值解实验报告

1、用有限元方法求下列边值问题的数值解：''()112x -y +y =2s i n ,0∈∈??∈(0,)?， 其中取1ν= 要求画出解曲面。迭代格式如下： 1221212111111111122142212n n n n n n j j j j j j n n n n n n j j j j j j V V V V V V h h V V V V V V h h τ++++++++++-+-??-()-()()-()??++?????? ??-+-+??=+??????

1、 %Ritz Galerkin方法求解方程 function u1=Ritz(x) %定义步长 h=1/100; x=0:h:1; n=1/h; a=zeros(n-1,1); b=zeros(n,1); c=zeros(n-1,1); d=zeros(n,1); %求解Ritz方法中内点系数矩阵 for i=1:1:n-1 b(i)=(1/h+h*pi*pi/12)*2; d(i)=h*pi*pi/2*sin(pi/2*(x(i)+h))/2+h*pi*pi/2*sin(pi/2*x(i+1))/2; end %右侧导数条件边界点的计算 b(n)=(1/h+h*pi*pi/12); d(n)=h*pi*pi/2*sin(pi/2*(x(i)+h))/2; for i=1:1:n-1 a(i)=-1/h+h*pi*pi/24; c(i)=-1/h+h*pi*pi/24; end %调用追赶法 u=yy(a,b,c,d) %得到数值解向量 u1=[0,u] %对分段区间做图 plot(x,u1) %得到解析解 y1=sin(pi/2*x); hold on plot(x,y1,'o') legend('数值解','解析解') function x=yy(a,b,c,d) n=length(b); q=zeros(n,1); p=zeros(n,1); q(1)=b(1); p(1)=d(1); for i=2:1:n

偏微分方程数值解

偏微分方程数值解 偏微分方程地构建科学、工程学和其他领域的数学模型的主要手段。一般情况下,这些模型都需要用数值方法去求解。本书提供了标准数值技术的简明介绍。借助抛物线型、双曲线型和椭圆型方程的一些简单例子介绍了常用的有限差分方法、有限元方法、有限体方法、修正方程分析、辛积分格式、对流扩散问题、多重网络、共轭梯度法。利用极大值原理、能量法和离散傅里叶分析清晰严格地处理了稳定性问题。本书全面讨论了这些方法的性质,并附有典型的图像结果,提供了不同难度的例子和练习。 本书可作为数学、工程学及计算机科学专业本科教材,也可供工程技术人员和应用工作者参考。 偏微分方程数值解---学习总结(2) 关于SobolveSobolve空间的几个重要定理 迹定理 : ΩΩ是 RdRd 的一个有界开子集，具有李普希茨连续边界?Ω?Ω, s>12s>12, 则 a.存在唯一的连续线性映射γ0:Hs(Ω)→Hs?12(?Ω),满足γ0v=v ∣∣?Ω,?v∈Hs(Ω)∩C0(Ωˉˉˉˉ), b.存在唯一的连续映射R0：Hs?12(?Ω)→Hs(Ω),满足γ0°R0°φ=φ,?φ∈Hs?12(?Ω).(1)(2)(1)a.存在唯一的连续线性映射γ0:Hs(Ω)→Hs?12(?Ω),满足γ0v=v|?Ω,?v∈

Hs(Ω)∩C0(Ωˉ),(2)b.存在唯一的连续映射R0：Hs?12(?Ω)→Hs(Ω),满足γ0°R0°φ=φ,?φ∈Hs?12(?Ω). 迹定理把区域内部与边界联系起来. 上面定理中边界?Ω?Ω当被它的一个子集ΣΣ代替时，结论依然成立. S=1时， γ0:H1(Ω)→H12(?Ω)?L2(?Ω)||γ0v||0,?Ω≤||γ0v||2,?Ω≤C||v||1=C(||v||0+||?v||0).γ0:H1(Ω)→H12(?Ω)? L2(?Ω)||γ0v||0,?Ω≤||γ0v||2,?Ω≤C||v||1=C(||v||0+||? v||0). 注意几个范数 ||?||k||?||0||?||1||??||0=||?||k,2=||?||L2=||?||1,2=(||?||20+||??||20)12=|?|1.(3)(4)(5)(6)(3)||?||k=||?||k,2(4)||? ||0=||?||L2(5)||?||1=||?||1,2=(||?||02+||??||02)12(6)||?? ||0=|?|1. 庞加莱不等式(Poincare inequality): 假设ΩΩ是 RdRd 的一个有界联通开子集，ΣΣ是边界?Ω?Ω的一个非空的李普希茨连续子集. 则存在一个常数 CΩ>0CΩ>0满足 ∫Ωv2(x)dx≤CΩ∫Ω|?v(x)|2dx,?v∈H1Σ(Ω),其中H1Σ(Ω)={v ∈H1(Ω),γΣv=v∣∣Σ=0}.∫Ωv2(x)dx≤CΩ∫Ω|?v(x)|2dx,?v∈HΣ1(Ω),其中HΣ1(Ω)={v∈H1(Ω),γΣv=v|Σ=0}.

常微分方程边值问题的数值解法

第8章 常微分方程边值问题的数值解法 引 言 第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。 只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见，我们以二阶边值问题为 则边值问题(8.1.1)有唯一解。 推论 若线性边值问题 ()()()()()(),, (),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤?? ==? (8.1.2) 满足 (1) (),()p x q x 和()f x 在[,]a b 上连续； (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。 求边值问题的近似解，有三类基本方法： (1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解； (2) 有限元法(finite element method)；

(3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。 差分法 8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法 设二阶线性常微分方程的边值问题为 (8.2.1)(8.2.2) ()()()(),,(),(), y x q x y x f x a x b y a y b αβ''-=<

偏微分方程数值解法试题与答案

一．填空（1553=?分） 1．若步长趋于零时，差分方程的截断误差0→lm R ，则差分方程的解lm U 趋近于微分方 程的解lm u . 此结论_______（错或对）； 2．一阶Sobolev 空间{} )(,,),()(21 Ω∈''=ΩL f f f y x f H y x 关于内积=1),( g f _____________________是Hilbert 空间； 3．对非线性（变系数）差分格式，常用 _______系数法讨论差分格式的_______稳定性； 4．写出3 x y =在区间]2,1[上的两个一阶广义导数：_________________________________， ________________________________________； 5．隐式差分格式关于初值是无条件稳定的. 此结论_______（错或对）。 二．（13分）设有椭圆型方程边值问题 用1.0=h 作正方形网格剖分 。 （1）用五点菱形差分格式将微分方程在内点离散化； （2）用截断误差为)(2 h O 的差分法将第三边界条件离散化； （3）整理后的差分方程组为 三．（12）给定初值问题 x u t u ??=?? ， ()10,+=x x u 取时间步长1.0=τ，空间步长2.0=h 。试合理选用一阶偏心差分格式（最简显格式）, 并以此格式求出解函数),(t x u 在2.0,2.0=-=t x 处的近似值。 1．所选用的差分格式是： 2．计算所求近似值： 四．（12分）试讨论差分方程 ()h a h a r u u r u u k l k l k l k l ττ + - = -+=++++11,111 1 逼近微分方程 0=??+??x u a t u 的截断误差阶R 。 思路一：将r 带入到原式，展开后可得格式是在点（l+1/2,k+1/2）展开的。 思路二：差分格式的用到的四个点刚好是矩形区域的四个顶点，可由此构造中心点的差分格 式。

偏微分方程数值解期末试题及答案(内容参考)

偏微分方程数值解试题(06B) 参考答案与评分标准 信息与计算科学专业 一（10分）、设矩阵A 对称，定义)(),(),(2 1 )(n R x x b x Ax x J ∈-= ，)()(0x x J λλ?+=.若0)0('=?,则称称0x 是)(x J 的驻点（或稳定点）.矩阵A 对称（不必正定），求证0x 是)(x J 的驻点的充要条件是：0x 是方程组 b Ax =的解 解: 设n R x ∈0是)(x J 的驻点,对于任意的n R x ∈,令 ),(2 ),()()()(2 000x Ax x b Ax x J x x J λλλλ?+ -+=+=, (3分) 0)0('=?,即对于任意的n R x ∈,0),(0=-x b Ax ,特别取b Ax x -=0,则有 0||||),(2000=-=--b Ax b Ax b Ax ,得到b Ax =0. (3分) 反之,若 n R x ∈0满足 b Ax =0,则对于任意的 x ,)(),(2 1 )0()1()(00x J x Ax x x J >+ ==+??,因此0x 是)(x J 的最小值点. (4分) 评分标准:)(λ?的展开式3分, 每问3分,推理逻辑性1分 二（10分）、 对于两点边值问题：????? ==∈=+-=0 )(,0)() ,()(' b u a u b a x f qu dx du p dx d Lu 其中]),([,0]),,([,0)(min )(]),,([0min ] ,[1b a H f q b a C q p x p x p b a C p b a x ∈≥∈>=≥∈∈ 建立与上述两点边值问题等价的变分问题的两种形式：求泛函极小的Ritz 形式和 Galerkin 形式的变分方程。 解: 设}0)(),,(|{11 =∈=a u b a H u u H E 为求解函数空间,检验函数空间.取),(1 b a H v E ∈,乘方程两端,积分应用分部积分得到 (3分) )().(),(v f fvdx dx quv dx dv dx du p v u a b a b a ==+=??,),(1 b a H v E ∈? 即变分问题的Galerkin 形式. (3分)

实验4常微分方程数值解

实验4 常微分方程数值解 化工系毕啸天2010011811 【实验目的】 1. 练习数值微分的计算； 2. 掌握用MATLAB 软件求微分方程初值问题数值解的方法； 3. 通过实例学习用微分方程模型解决简化的实际问题； 4. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。 【实验内容】 题目3 小型火箭初始重量为1400kg，其中包括1080kg 燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N 的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。 3.1 燃料燃烧过程物理模型分析 设火箭质量为m，高度为h，速度为v，加速度为a，火箭推力为F，重力加速度为g，阻力为f。 1.由火箭上总共携带燃料1080kg，燃料燃烧率为18kg/s，可知火箭上升时间t=60s时，燃料全部烧尽。 2.由阻力正比于速度的平方，比例系数0.4kg/m，可知阻力表达式为f=0.4v2。 3.由于燃料燃烧，火箭的质量是时间的函数，易知m(t)=m0-18t 4. 5.根据牛顿第二运动定律，有。代入数据有 解出 由以上5条分析，我们得到了一个常微分方程组： 初值条件为：v0=0,h(0)=0,t60s. 3.2 程序代码 根据常微分方程组的初值问题，在MATLAB中计算数值解。 记x(1) = h，x(2)= v，x =(x(1), x(2))T 首先编写M文件 function dx = Rocket(t,x) dx=[x(2);(32000-0.4*x(2)^2)/(1400-18*t)-9.8];%以向量形式表示微分方程 end ts=0:60 %终点时间为60s，步长定义为1即可

偏微分方程数值解法

一、 问题 用有限元方法求下面方程的数值解 2 u u u f t ?-?+=? in (]0,T Ω? 0u = on []0,T ?Ω? ()00,u x u = in Ω 二、 问题分析 第一步 利用Green 公式，求出方程的变分形式 变分形式为：求()()21 00,;u L T H ∈Ω，使得 ()())(2 ,,,,u v u v u v f v t ???+??+= ???? ()10v H ?∈Ω （*） 以及 ()00,u x u =. 第二步 对空间进行离散，得出半离散格式 对区域Ω进行剖分，构造节点基函数，得出有限元子空间：()12,,,h NG V span ???=???，则（*）的Galerkin 逼近为： []0,t T ?∈，求()()1 0,h h u t x V H ∈?Ω，使得 ()()()()() () )(2 ,,,,h h h h h h h d u t v u t v u t v f v dt +??+= h h v V ?∈ （**） 以及()0,0h h u u =，0,h u 为初始条件0u 在h V 中的逼近，设0,h u 为0u 在h V 中的插值. 则0t ?≥，有()()1 N G h i i i u t t ξ? == ∑，0,h u =01 N G i i i ξ?=∑，代人（**）即可得到一常微分方程组. 第三步 进一步对时间进行离散，得到全离散的逼近格式 对 du dt 用差分格式.为此把[]0,T 等分为n 个小区间[]1,i i t t -，其长度1i i T t t t n -?=-= ,n t T =. 这样把求i t 时刻的近似记为i h u ，0 h u 是0u 的近似.这里对（**）采用向后的欧拉格式，即 ()()() () )(2 11 11 1 ,,,,i i i i h h h h h h h i h u u v u v u v f v t ++++-+??+ = ? h h v V ?∈ (***) i=0,1,2…,n-1. 0 h u =0,h u 由于向后欧拉格式为隐式格式且含有非线性项，故相邻两时间步之间采用牛顿迭代，即：

偏微分方程数值解复习题(2011硕士)

偏微分方程数值解期末复习（2011硕士） 一、考题类型 本次试卷共六道题目，题型及其所占比例分别为： 填空题20%；计算题80% 二、按章节复习内容 第一章 知识点：Euler法、向前差商、向后差商、中心差商、局部截断误差、整体截断误差、相容性、收敛性、阶、稳定性、显格式、隐格式、线性多步法、第一特征多项式、第二特征多项式、稳定多项式、绝对稳定等； 要求: 会辨认差分格式, 判断线性多步法的误差和阶; 第二章 知识点：矩形网格、(正则,非正则)内点、边界点、偏向前(向后,中心)差商、五点差分格式、增设虚点法、积分插值法、线性椭圆型差分格式、极值原理、比较定理、五点差分格式的相容收敛和、稳定性等； 要求: 建立椭圆型方程边值问题的差分格式, 极值原理; 第四章 知识点：最简显格式、最简隐格式、CN格式、双层加权格式、Richardson 格式、网格比、传播因子法(分离变量法) 、传播因子、传播矩阵、谱半径、von Neumann条件、跳点格式、ADI格式、线性椭圆型差分格式、极值原理、比较定理、五点差分格式的相容收敛和稳定性等； 要求: 建立抛物型方程边值问题的差分格式, 计算局部截断误差; 第五章 知识点：左偏心格式、右偏心格式、中心格式、LF格式、LW格式、Wendroff 格式、跳蛙格式、特征线、CFL条件等； 要求: 建立双曲型方程边值问题的差分格式, 计算局部截断误差; 第七章 要求: 会用线性元(线性基)建立常微分方程边值问题的有限元格式

三 练习题 1、 已知显格式21131()22 n n n n u u h f f +++-=-，试证明格式是相容的，并求它的阶。 P39+P41 2、用Taylor 展开原理构造一元函数一阶导数和二阶导数的数值微分公式。 提示：向前、向后和中心差商与一阶导数间关系，二阶中心差商与二阶导数 之间的关系 课件 3、用数值微分方法或数值积分方法建立椭圆型方程 2222(,),(,),u u f x y x y x y ??--=?∈Ω?? :01,01x y Ω≤≤≤≤ 内点差分格式。 P75+课件 4、构造椭圆型方程边值问题的差分格式. P101 (4)题 5、构建一维热传导方程220,(0)u u Lu a a t x ??=-=>??的数值差分格式(显隐格式等)。 参考P132-135相关知识点 6、设有逼近热传导方程22(0)u u Lu a f a const t x ??≡-==>??的带权双层格式 ()()1111111122(1)2k k j j k k k k k k j j j j j j u u a u u u u u u h θθτ++++-+-+-??=-++--+?? 其中[0,1]θ∈，试求其截断误差。并证明当2 1212h a θτ=-时，截断误差的阶最 高阶为24()O h τ+。 P135+P165+课件 7、传播因子法证明抛物型方程22(0)u u Lu a f a const t x ??≡-==>??的最简显隐和六点CN 格式稳定性。 P156+课件 8、对一阶常系数双曲型方程的初边值问题 0,0,0,0,(,0)(),0,(0,)(),0, u u a t T x a t x u x x x u t t t T φψ???+=<≤<<∞>?????=≤<∞??=≤≤?

常微分方程数值解法

i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。

偏微分方程数值解法

“十二五”国家重点图书出版规划项目 信息与计算科学丛书 67 偏微分方程数值解法 陈艳萍鲁祖亮刘利斌编著

内 容 简 介 本书试图用较少的篇幅描述偏微分方程的几种数值方法. 主要内容包括：Sobolev空间初步, 椭圆边值问题的变分问题, 椭圆问题的有限差分方法, 抛物型方程的有限差分方法, 双曲型方程的有限差分方法, 椭圆型方程的有限元方法, 抛物及双曲方程的有限元方法, 椭圆型方程的混合有限元方法, 谱方法等. 本书内容丰富, 深入浅出, 尽可能地用简单的方法来描述一些理论结果, 并根据作者对有限差分、有限元、混合有限元、谱方法的理解和研究生教学要求, 全面、客观地评价各种数值计算方法，并列举一些数值计算的例子, 阐述许多新的学术观点. 本书可作为高等学校数学系高年级本科生和研究生的教材或参考书, 也可作为计算数学工作者和从事科学与工程计算的科研人员的参考书. 图书在版编目(CIP)数据 偏微分方程数值解法/陈艳萍, 鲁祖亮, 刘利斌编著. —北京：科学出版社, 2015.1 (信息与计算科学丛书67) ISBN 978-7-03-000000-0 Ⅰ. ①偏… Ⅱ. ①陈… ②鲁… ③刘… Ⅲ. ① Ⅳ.① 中国版本图书馆CIP数据核字(2014) 第000000号 责任编辑: 王丽平／责任校对: 彭涛 责任印制: 肖钦／封面设计: 陈敬 出版 北京东黄城根北街16号 邮政编码: 100717 https://www.wendangku.net/doc/4b16914727.html, 印刷 科学出版社发行 各地新华书店经销 * 2015年1月第一版开本: 720×1000 1/16 2015年1月第一次印刷印张: 14 字数: 280 000 定价: 88.00元 (如有印装质量问题, 我社负责调换)

matlab常微分方程的数值解法实验报告

实验四 常微分方程的数值解法 指令： [t,y]=ode23(‘fun ’,tspan,yo) 2/3阶龙格库塔方法 [t,y]=ode45(‘fun ’,tspan,yo) 4/5阶龙格库塔方法 [t,y]=ode113(‘fun ’,tspan,yo) 高阶微分方程数值方法 其中fun 是定义函数的文件名。该函数fun 必须以为dx 输出量，以t,y 为输入量。tspan=[t0 tfina]表示积分的起始值和终止值。yo 是初始状态列向量。 考虑到初始条件有 00 d , (0)0,d d , (0)0. d S SI S S t I SI I I I t ββμ?=-=>??? ?=-=≥?? (5.24) 这就是Kermack 与McKendrick 的SIR 仓室模型. 方程(5.24)无法求出()S t 和()I t 的解析解.我们先做数值计算。Matlab 代码为： function dy=rigid(t,y) dy=zeros(2,1); a=1; b=0.3; dy(1)=a*y(1).*y(2)-b*y(1); dy(2)=-a*y(1).*y(2); ts=0:.5:50; x0=[0.02,0.98]; [T,Y]=ode45('rigid',ts,x0); %plot(T,Y(:,1),'-',T,Y(:,2),'*') plot(Y(:,2),Y(:,1),'b--') xlabel('s') ylabel('i') 任务： 1 画出i （t ）， 2分析各参数的影响 例57：求解两点边值问题：0)5(,0)1(,32 ==='-''y y x y y x 。（注意：相应的数值解法比较复杂）。 y=dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=0','x') ↙ y = -1/3*x^3+125/468+31/468*x^4

偏微分方程数值解法答案

1. 课本2p 有证明 2. 课本812,p p 有说明 3. 课本1520,p p 有说明 4. Rit2法，设n u 是u 的n 维子空间，12,...n ???是n u 的一组基底，n u 中的任一元素n u 可 表为1n n i i i u c ?==∑ ，则,11 11()(,)(,)(,)(,)22j n n n n n n i j i j j i j j J u a u u f u a c c c f ???=== -=-∑∑是12,...n c c c 的二次函数，(,)(,)i j j i a a ????=，令 () 0n j J u c ?=?，从而得到12,...n c c c 满足1 (,)(,),1,2...n i j i j i a c f j n ???===∑，通过解线性方程组，求的i c ，代入1 n n i i i u c ?==∑， 从而得到近似解n u 的过程称为Rit2法 简而言之，Rit2法：为得到偏微分方程的有穷维解，构造了一个近似解，1 n n i i i u c ?== ∑， 利用,11 11()(,)(,)(,)(,)22j n n n n n n i j i j j i j j J u a u u f u a c c c f ???===-=-∑∑确定i c ，求得近似解n u 的过程 Galerkin 法：为求得1 n n i i i u c ? == ∑形式的近似解，在系数i c 使n u 关于n V u ∈，满足(,)(,) n a u V f V =，对任 意 n V u ∈或（取 ,1j V j n ?=≤≤） 1 (,)(,),1,2...n i j i j i a c f j n ???===∑的情况下确定i c ，从而得到近似解1 n n i i i u c ?==∑的过程称 Galerkin 法为 Rit2-Galerkin 法方程： 1 (,)(,)n i j i j i a c f ???==∑ 5. 有限元法：将偏微分方程转化为变分形式，选定单元的形状，对求解域作剖分，进而构 造基函数或单元形状函数，形成有限元空间，将偏微分方程转化成了有限元方程，利用 有效的有限元方程的解法，给出偏微分方程近似解的过程称为有限元法。 6. 解：对求解区间进行网格剖分，节点01......i n a x x x x b =<<<<=得到相邻节点1,i i x x -

微分方程数值解法实验报告

微分方程数值解法实验报告 姓名： 班级： 学号：

一：问题描述 求解边值问题： ()2(sin cos cos sin (0,1)(0,1)0,(,)x y u e x y x y G u x y G ππππππ+???=+??∈=?????=∈??? （x,y) 其精确解为)sin()sin(),()(y x e y x u y x πππ+= 问题一：取步长h=k=1/64,1/128,作五点差分格式，用Jacobi 迭代法，Gauss_Seidel 迭代法，SOR 迭代法（w=1.45）。求解差分方程，以前后两次重合到小数点后四位的迭代值作为解的近似值，比较三种解法的迭代次数以及差分解)128/1,64/1)(,(=h y x u h 与精确解的精度。 问题二：取步长h=k=1/64,1/128,作五点差分格式，用单参数和双参数PR 法解差分方程，近似到小数点后四位。与SOR 法比较精度和迭代步数。 问题三：取步长h=k=1/64,1/128,作五点差分格式，用共轭梯度法和预处理共轭梯度法解差分方程，近似到小数点后四位。与SOR 法与PR 法比较精度和迭代步数。 二．实验目的： 分别使用五点差分法（Jacobi 迭代，Gauss_Seidel 迭代，SOR 迭代），PR 交替隐式差分法（单参数，双参数），共轭梯度法，预共轭梯度法分别求椭圆方程的数值解。 三．实验原理： (1)Jacobi 迭代法 设线性方程组 (1) 的系数矩阵A 可逆且主对角元素均不为零,令 b Ax =nn a ,...,a ,a 2211

并将A 分解成 (2) 从而(1)可写成 令 其中. (3) 以为迭代矩阵的迭代法(公式) (4) 称为雅可比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为 (5) 其中 为初始向量. (2) Guass-Seidel 迭代法 由雅可比迭代公式可知,在迭代的每一步计算过程中是用的 全部分量来计算的所有分量,显然在计算第i 个分量时,已 经计算出的最新分量没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第 次近似的分量加以利用,就得到所谓解方程组的高斯 —塞德（Gauss-Seidel ）迭代法. 把矩阵A 分解成 (6) 其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成 即 其中 (7) 以为迭代矩阵构成的迭代法(公式) (8) ()nn a ,...,a ,a diag D 2211=()D D A A +-=() b x A D Dx +-=11f x B x +=b D f ,A D I B 1111 --=-=1B ()()111f x B x k k +=+? ??[],...,,k ,n ,...,i x a b a x n i j j )k (j j i i ii )k (i 21021111==∑-=≠=+()()()()()T n x ,...x ,x x 002010=()k x ()1+k x ()1+k i x ()() 1111+-+k i k x ,...,x 1+k ()1+k x ()1+k j x U L D A --=()nn a ,...,a ,a diag D 2211=U ,L --A ()b Ux x L D +=-22f x B x +=()()b L D f ,U L D B 1212---=-=2B ()()221f x B x k k +=+

第十章-偏微分方程数值解法

第十章 偏微分方程数值解法 偏微分方程问题，其求解十分困难。除少数特殊情况外，绝 大多数情况均难以求出精确解。因此，近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。 §1 差分方法的基本概念 1.1 几类偏微分方程的定解问题 椭圆型方程：其最典型、最简单的形式是泊松（Poisson ）方程 ),(22 2 2y x f y u x u u =??+??=? 特别地，当 0),(≡y x f 时，即为拉普拉斯（Laplace ）方程，又称 为调和方程 22 22 =??+??=?y u x u u Poisson 方程的第一边值问题为 ?? ?? ?Ω ?=Γ=Ω∈=??+??Γ∈),(),(),(),(),(22 22y x y x u y x y x f y u x u y x ? 其中 Ω为以Γ为边界的有界区域，Γ为分段光滑曲线， ΓΩY 称为定解区域，),(y x f ，),(y x ?分别为Ω，Γ上的已知连 续函数。 第二类和第三类边界条件可统一表示为

),(),(y x u u y x ?α=??? ? ??+??Γ∈n 其中n 为边界Γ的外法线方向。当0=α时为第二类边界条件， 0≠α时为第三类边界条件。 抛物型方程：其最简单的形式为一维热传导方程 2 20(0)u u a a t x ??-=>?? 方程可以有两种不同类型的定解问题： 初值问题 ?? ???+∞ <<∞-=+∞<<-∞>=??-??x x x u x t x u a t u )()0,(,00 22 ? 初边值问题 2 212 00,0(,0)()0(0,)(),(,)()0u u a t T x l t x u x x x l u t g t u l t g t t T ????-=<<<

常微分方程初值问题数值解法

常微分方程初值问题数值解法 朱欲辉 （浙江海洋学院数理信息学院, 浙江舟山316004） [摘要]：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法.然而在生产实 际和科学研究中所遇到的问题往往很复杂, 在很多情况下都不可能给出解的解析表达式. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler法、改进的Euler法、Runge－Kutta法以及线性多步法中的Adams显隐式公式和预测校正 公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 分别用不同的方法计算出近似解, 从得出的结果对比各种方法的优缺点. [关键词]：常微分方程;初值问题; 数值方法; 收敛性; 稳定性; 误差估计 Numerical Method for Initial-Value Problems Zhu Yuhui (School of Mathematics, Physics, and Information Science, Zhejiang Ocean University, Zhoushan, Zhejiang 316004) [Abstract]:In the course about ordinary differential equations, the methods for analytic solutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can’t be e xpressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis. [Keywords]:Ordinary differential equation; Initial-value problem; Numerical method; Convergence; Stability;Error estimate

常微分方程数值解法

第八章 常微分方程数值解法 考核知识点： 欧拉法，改进欧拉法，龙格-库塔法，单步法的收敛性与稳定性。 考核要求： 1． 解欧拉法，改进欧拉法的基本思想；熟练掌握用欧拉法，改进欧拉法、求微 分方程近似解的方法。 2． 了解龙格-库塔法的基本思想；掌握用龙格-库塔法求微分方程近似解的方 法。 3． 了解单步法的收敛性、稳定性与绝对稳定性。 例1 用欧拉法，预估——校正法求一阶微分方程初值问题 ? ??=-='1)0(y y x y ，在0=x （0，1）0.2近似解 解 （1）用1.0=h 欧拉法计算公式 n n n n n n x y y x y y 1.09.0)(1.01+=-+=+，1.0=n 计算得 9.01=y 82.01.01.09.09.02=?+?=y （2）用预估——校正法计算公式 1,0)(05.01.09.0)0(111)0(1=???-+-+=+=++++n y x y x y y x y y n n n n n n n n n 计算得 91.01=y ，83805.02=y 例2 已知一阶初值问题 ???=-='1 )0(5y y y 求使欧拉法绝对稳定的步长n 值。 解 由欧拉法公式 n n n n y h y h y y )51(51-=-=+ n n y h y ~)51(~1-=+

相减得01)51()51(e h e h e n n n -==-=-Λ 当 151≤-h 时，4.00≤

偏微分方程的数值解法

《偏微分方程数值解法》试题 (专业:凝聚态物理学号:2013201260 姓名:鄢建军) 1.考虑定解问题 （1）用迎风格式(P、45)求解 1,0 (,0) 0,0 t x u u x u x x += ? ? ≤ ? ? =? ?> ? ? 。 利用迎风格式编写Fortran程序语言,运行结果如下: Fig 1、迎风格式求解结果 （2）用Beam-Warming格式(P、51)求解。 利用Beam—Warming格式编写Fortran程序语言,运行结果如下 :

Fig 2、 Beam —Warming 格式求解结果 （3） 比较两种方法结果的异同。 将两种格式运行的结果绘制在一起,要求时间步长与空间步长在两种格式中都相同,运行结果如下图所示: Fig 3、 迎风格式与Beam-Warming 格式求解结果比较 从两种格式的运行结果来瞧,都存在边缘的误差现象,相比而言,Beam-Warming 格式的运行结果差一些。但就是理论上分析,迎风格式的截断误差为()h οτ+,而Beam-Warming 格式的截断误差为22()h h οττ++。稳定性上来分析,迎风格式的稳定性较好,要求1(/)a h λλτ≤=,Beam-Warming 格式的稳定性条件为2(/)a h λλτ≤=。 2. 考虑定解问题212 1110,04(,0)sin ,0(0,)(,)0u u a x l t t u x x x l l u t u l t π???-=<