第4章常微分方程数值解法PPT课件

合集下载

常微分方程的数值解法

常微分方程的数值解法

常微分方程的数值解法在自然科学的许多领域中,都会遇到常微分方程的求解问题。

然而,我们知道,只有少数十分简单的微分方程能够用初等方法求得它们的解,多数情形只能利用近似方法求解。

在常微分方程课中已经讲过的级数解法,逐步逼近法等就是近似解法。

这些方法可以给出解的近似表达式,通常称为近似解析方法。

还有一类近似方法称为数值方法,它可以给出解在一些离散点上的近似值。

利用计算机解微分方程主要使用数值方法。

我们考虑一阶常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(yx y y x f dx dy在区间[a, b]上的解,其中f (x, y )为x, y 的已知函数,y 0为给定的初始值,将上述问题的精确解记为y (x )。

数值方法的基本思想是:在解的存在区间上取n + 1个节点b x x x x a n =<<<<= 210这里差i i i x x h -=+1,i = 0,1, …, n 称为由x i 到x i +1的步长。

这些h i 可以不相等,但一般取成相等的,这时na b h -=。

在这些节点上采用离散化方法,(通常用数值积分、微分。

泰勒展开等)将上述初值问题化成关于离散变量的相应问题。

把这个相应问题的解y n 作为y (x n )的近似值。

这样求得的y n 就是上述初值问题在节点x n 上的数值解。

一般说来,不同的离散化导致不同的方法。

§1 欧拉法与改进欧拉法 1.欧拉法1.对常微分方程初始问题(9.2))((9.1) ),(00⎪⎩⎪⎨⎧==y x y y x f dx dy用数值方法求解时,我们总是认为(9.1)、(9.2)的解存在且唯一。

欧拉法是解初值问题的最简单的数值方法。

从(9.2)式由于y (x 0) = y 0已给定,因而可以算出),()('000y x f x y =设x 1 = h 充分小,则近似地有:),()(')()(00001y x f x y hx y x y =≈-(9.3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为y (x 1)的近似值。

常微分方程的数值解

常微分方程的数值解

f ( x, y1 ) f ( x, y2 ) L y1 y2
(其中 L 为 Lipschitz 常数)则初值问题( 1 )存 在唯一的连续解。
求问题(1)的数值解,就是要寻找解函数在一 系列离散节点x1 < x2 <……< xn < xn+1 上的近似 值y1, y 2,…,yn 。 为了计算方便,可取 xn=x0+nh,(n=0,1,2,…), h称为步长。
(1),(2)式称为初值问题,(3)式称为边值问题。 在实际应用中还经常需要求解常微分方程组:
f1 ( x, y1 , y2 ) y1 ( x0 ) y10 y1 (4) f 2 ( x, y1 , y2 ) y2 ( x0 ) y20 y2
本章主要研究问题(1)的数值解法,对(2)~(4)只 作简单介绍。
得 yn1 yn hf ( xn1 , yn1 )
上式称后退的Euler方法,又称隐式Euler方法。 可用迭代法求解
二、梯形方法 由
y( xn1 ) y( xn )
xn1 xn
f ( x, y( x))dx
利用梯形求积公式: x h x f ( x, y( x))dx 2 f ( xn , y( xn )) f ( xn1 , y( xn1 ))
常微分方程的数言 简单的数值方法 Runge-Kutta方法 一阶常微分方程组和高阶方程
引言
在高等数学中我们见过以下常微分方程:
y f ( x, y, y) a x b y f ( x, y ) a x b (2) (1) (1) y ( x ) y , y ( x ) y 0 0 0 0 y ( x0 ) y0 y f ( x, y, y) a x b (3) y(a) y0 , y(b) yn

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

常微分方程初值问题数值解法
根据微分方程的性质和初始条件,常 微分方程初值问题可以分为多种类型, 如一阶、高阶、线性、非线性等。
数值解法的必要性
实际应用需求
许多实际问题需要求解常微分方程初值问题,如物理、 化学、生物、工程等领域。
解析解的局限性
对于复杂问题,解析解难以求得或不存在,因此需要 采用数值方法近似求解。
数值解法的优势
未来发展的方向与挑战
高精度算法
研究和发展更高精度的算法,以提高数值解的准确性和稳定性。
并行计算
利用并行计算技术,提高计算效率,处理大规模问题。
自适应方法
研究自适应算法,根据问题特性自动调整计算精度和步长。
计算机技术的发展对数值解法的影响
1 2
硬件升级
计算机硬件的升级为数值解法提供了更强大的计 算能力。
它首先使用预估方法(如欧拉方法)得到一个 初步解,然后使用校正方法(如龙格-库塔方法) 对初步解进行修正,以提高精度。
预估校正方法的优点是精度较高,且计算量相 对较小,适用于各种复杂问题。
步长与误差控制
01
在离散化过程中,步长是一个重要的参数,它决定 了离散化的精度和计算量。
02
误差控制是数值逼近的一个重要环节,它通过设定 误差阈值来控制计算的精度和稳定性。
能够给出近似解的近似值,方便快捷,适用范围广。
数值解法的历史与发展
早期发展
早在17世纪,科学家就开始尝 试用数值方法求解常微分方程。
重要进展
随着计算机技术的发展,数值 解法在20世纪取得了重要进展, 如欧拉法、龙格-库塔法等。
当前研究热点
目前,常微分方程初值问题的 数值解法仍有许多研究热点和 挑战,如高精度算法、并行计
软件优化
软件技术的发展为数值解法提供了更多的优化手 段和工具。

常微分方程数值解法5262115页PPT文档

常微分方程数值解法5262115页PPT文档
x 1 ( t ) 表示时刻 t 食饵的密度,x 2 ( t ) 表示捕食者的密度;
r 表示食饵独立生存时的增长率;
d 表示捕食者独立生存时的死亡率;
a 表示捕食者的存在对食饵增长的影响系数,反映捕
食者对食饵的捕获能力;
b 表示食饵的存在对捕食者增长的促进系数,反映食
饵对捕食者的喂养能力
150 100
令 y 1 y ,y 2 y ',y 3 y '', ,y n y ( n 1 )
可以将以上高阶微分方程化为如下一阶常微分方程组
y1 ' y2 y2 ' y3 yn ' an(x)y1
a1(x)yn f (x)
例:P120,1(a),Bessel方程
常微分方程的数值解
一般地,凡表示未知函数,未知函数的导 数与自变量之间的关系的方程叫做微分方 程.未知函数是一元函数的,叫常微分方 程;未知函数是多元函数的,叫做偏微分方 程.

y ' x y'x2y2 y''y'xy
Matlab实现 [t,x]=ode45(f,ts,x0,options,p1,p2,......)
50 0 0
30 20 10
0 0
10
20
50
30
20
10

0
30
0
10
8
6
4
2
100
0
50
100
150
50
100
高阶常微分方程的解法
高阶常微分方程
y ( n ) a 1 ( x ) y ( n 1 ) a ( n 1 ) ( x ) y ' a n ( x ) y f( x )

常微分方程 PPT课件

常微分方程 PPT课件

分曲线和积分曲线族的概念,只不过此时积分曲线所在的空间维数不同,我们将
在第4章详细讨论.
最后,我们要指出,本书中按习惯用
代替
而 分别代表
本节要点: 1.常微分程的定义,方程的阶,隐式方程,显式方程,线性方程,非线性方程. 2.常微分方程解的定义,通解,特解,通积分,特积分. 3.初值问题及初值问题解的求法. 4.解的几何意义,积分曲线.
所以它们都是一阶齐次方程.因此,一阶齐次微分方程可以 写为
(1.27)
1.3.1 齐次方程的解法 方程(1.27)的特点是它的右端是一个以为
变元的函数,经过如下的变量变换,它能化 为变量可分离方程.
令 则有 代入方程(1.27)得
(1.28)

方程(1.28)是一个 变量可分离方程,当 时,分离 变量并积分,得到它的通积分 (1.29)
常微分方程课件
第一章 初等积方法 第二章 基本定理 第三章 线性微分方程 第四章 线性微分方程组 第五章 定性与稳定性概念 第六章 一阶偏微方程初步
第1讲 微分方程与解 微分方程
什么是微分方程?它是怎样产生的?这是首先要回答的问题.
300多年前,由牛顿(Newton,1642-1727)和莱布尼兹 (Leibniz,1646-1716)所创立的微积分学,是人类科学史上划
(1.13)
显然,方程(1.4)是一阶线性方程;方程(1.5)是一阶非线性方程;方程 (1.6)是二阶线性方程;方程(1.7)是二阶非线性方程.
通解与特解
微分方程的解就是满足方程的函数,可定义如下.
定义1.1 设函数 在区间I上连续,且有直
到n阶的导数.如果把
代入方程(1.11),得到在
区间I上关于x的恒等式,

实验4 常微分方程的数值解法

实验4 常微分方程的数值解法

[内容] 1. 欧拉格式(或改进的欧拉格式),编写 相应的程序并能正确运行。 2. 经典四:先描述清楚问题。
实验4 常微分方程的数值解法 [要求]
1.程序的调试要耐心、细致;
2.语句应尽可能加注注释; 3.本次实验的各个程序(M文件)打包成压缩文件 (格式:学号姓名.RAR,如:200910119李娟.RAR), 按时提交。
实验4 常微分方程的数值解法
[目的]
1.常微分方程差分算法的计算机实现;
2.进一步理解欧拉格式、改进的欧拉格式(预报—
校正系统)、龙格—库塔格式等算法,会运用这些方
法解决初步的常微分方程的求解问题; 3.进一步熟悉MATLAB数学软件的使用,锻炼程序 调试、排错的能力。
实验4 常微分方程的数值解法

常微分方程数值解法

常微分方程数值解法

欧拉方法
总结词
欧拉方法是常微分方程数值解法中最基础的方法之一,其基本思想是通过离散化时间点上的函数值来 逼近微分方程的解。
详细描述
欧拉方法基于微分方程的局部线性化,通过在时间点上逐步逼近微分方程的解,得到一系列离散点上 的近似值。该方法简单易行,但精度较低,适用于求解初值问题。
龙格-库塔方法
总结词
影响
数值解法的稳定性对计算结果的精度和可靠 性有重要影响。
判断方法
通过分析数值解法的迭代公式或离散化方法, 判断其是否具有稳定性和收敛性。
数值解法的收敛性
定义
数值解法的收敛性是指随着迭代次数的增加, 数值解逐渐接近于真实解的性质。
影响
数值解法的收敛性决定了计算结果的精度和 计算效率。
分类
根据收敛速度的快慢,可以分为线性收敛和 超线性收敛等。
判断方法
通过分析数值解法的迭代公式或离散化方法, 判断其是否具有收敛性。
误差分析
定义
误差分析是指对数值解法计算过程中 产生的误差进行定量分析和估计的过 程。
分类
误差可以分为舍入误差、截断误差和 初始误差等。
影响
误差分析对于提高计算精度和改进数 值解法具有重要意义。
分析方法
通过建立误差传递公式或误差估计公 式,对误差进行定量分析和估计。
生物学
生态学、生物种群动态和流行病传播 等问题可以通过常微分方程进行建模
和求解。
化学工程
化学反应动力学、化学工程流程模拟 等领域的问题可以通过常微分方程进 行描述和求解。
经济学
经济系统动态、金融市场模拟和预测 等问题可以通过常微分方程进行建模 和求解。
02 常微分方程的基本概念
常微分方程的定义

常微分方程的数值解及其它问题PPT课件

常微分方程的数值解及其它问题PPT课件

2020/10/13
3
• 数值积分
• quad('fname',a,b,tol,trace) Simpson法求数值积 分。
• quad8('fname',a,b,tol,trace) Newton-Cotes法求 数值积分。
• fname是被积函数文件名。
• b,a分别是积分上下限。
• 用tol来控制积分精度。可缺省。缺省时默认 tol=0.001。
b
1
( b a )n
a f ( x ) d ( b x a ) 0 f ( a ( b a ) u ) d u ni 1 f ( a ( b a ) u i)
2020/10/13
13
例 如 对 于 上 面 的 离 散 函 数 , 我 们 有 a=1 , b=1.5。于是我们可键入:
• 用trace来控制是否用显示积分过程。可缺省。 缺省时默认trace=0,不显示过程。
2020/10/13
4
例如:求
e3 x2
0
dx。
第一步:在编辑器中建立被积函数的M文
件。取名为fname。即在编辑器中输入:
functiห้องสมุดไป่ตู้n y=fname(x)
y=exp(-x^2);
2020/10/13
例如有函数
x i 1 . 0 0 0 0 1 . 1 0 0 0 1 . 2 0 0 0 1 . 3 0 0 0 1 . 4 0 0 0 1 . 5 0 0 0 y i 1 . 8 4 1 5 2 . 1 9 0 3 2 . 5 5 8 4 2 . 9 4 2 6 3 . 3 3 9 6 3 . 7 4 6 2
常微分方程的数值解及其它问 题的数值解

微分方程数值解

微分方程数值解

第四章 微分方程数值解4.1 算 法当常微分方程能解析求解时,可利用Matlab 符号工具箱中的功能找到精确解. 见下例求解方程20y y y '''++=. 键入:syms x y %定义符号变量 diff_equ= ‘D2y+2*Dy -y=0’; %D2y 表示,y ''Dy=y 'y=dsolve (diff_equ, ‘x’)%定义x 为自变量y=cl*exp ((2^(1/2)-1)*x+c2*exp (-(2^(1/2)+1)*x) %表达式中含c1与c2,表示通解.%初始条件为y (0)=0,y '(0)=1时,按如下方式调用y=dsolve (diff_equ, ‘y (0)=0’, ‘Dy (0)=1’, ‘x’)y=1/4*2^(1/2)*exp ((2^(1/2)-1)*x)— 1/4*2^(1/2)*exp (-(2^(1/2)+1)*x) %画出函数y=y (x)的图形ezplot (y,[-2,2])图形具体形式请上机试之.在方程无法获得解析解的情况下,可方便地获得数值解. 下面的例子说明用Matlab 求数值解的方法及应注意的问题.[例1] 求解范德堡(vander pol )方程222(1)0d x dx x x dt dtμ--+=求解高阶方程,必须等价地变换为一阶微分方程组,对本例,通过定义两个新的变量,实现这一变换1,2/y x y dx dt ==则令 1/2dy dt y =22/(11)*21dy dt y y y μ=--编写求解程序分为两部分,第一部分为待求解的方程,存盘的文件名为<待求解方程的函数名.m >,第二部分为求解主程序,本例中取名为main1.m.首先编写待求解方程的M -文件. 文件存盘名为“vdpol.m ”. function yprime=vdpol (,)t y yprime (1)=y (2);mu=2yprime (2)=mu*(1(1)^2)*(2)(1)y y y --;yprime=[yprime (1);yprime (2)];说明 函数yprime=vdpol (,)t y 中. 定义t 为自变量,y 的形式取决于求解方程的阶数,本例中,[(1),(2)],(1)y y y y =为解向量,(2)y 为导数向量. yprime (1)(1)y '=, yprime (2)(1)y ''=,函数返回vander pol 方程的导数列向量. 因为所求结果为方程数值解,所以各向量维数只有在主程序求解时定下精度后才能确定.主程序定名为main1.m ,你可用你所喜欢的其它名子,但vdpol.m 除外. clear functions%调试程序时,放置这一语句是必要的. 它清除前边已编译的存在于内存中的废弃程序 [,t y ]=ode23 (‘vdpol’,[0,30],[1,0]); y1=y (:,1); %解曲线. y2=y (:,2); %解曲线的导数. polt (,1,,2,t y t y ‘_ _’)说明 龙格_库塔的2阶与4阶改进型求解公式的实现,其指令分别为: [,t x ]=ode23 (‘f ’,,0,ts x options) [,t x ]=ode45 (‘f ’,,0,ts x options)其中ts 可由系统依据精度要求自动设定,亦可由使用者依据实际需要自己确定,分别说明之. (1)若令[0,1,,]ts t t tf =,则输出在指定时刻0,1,,t t tf 给出,当0::ts t k tf =时,输出在区间[0,]t tf 的等分点上给出,k 为步长.(2)若[0,],0ts t tf t =为自变量初值,tf 为终值,此时,options 决定自变量t 的维数,t 中的时间点不是等间隔的,这是为了保证所需的相对精度,积分算法改变了步长. 用于设定误差限的参数options 可缺省,此时系统设定相对误差为310-,绝对误差为610-,若自行设定误差限,可用如下语句: options=odeset (‘reltol’,rt , ‘abstol’,at ) 这里的rt 与at 分别为设定的相对与绝对误差.须注意的是无论用哪种方法确定t 的取值方式,0,t tf 必须由使用者确定且应与0x 相匹配. 0x 为初始条件,本例中0[1,0]y =,因为[0,30]ts =,这意味着解曲线(0)1y =,(0)0y '=,一般说,当解n 个未知函数的方程组时,0x 为n 维向量,共含有n 个初始条件.两个输出参数是列向量t 与矩阵x ,它们具有相同的行数,而矩阵x 的列数等于方程组的个数,本例中y 的列数为2,其中,(:,1)y 为自变量t 上各点函数值,(:,2)y 为t 上各点导数值.最后,提请读者注意的是:ode45也不总是比ode23好,在很多时候,低阶算法更有效,有关微分方程数值解法的更进一步信息,请参考数值分析方面的书籍. 有些参考书提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法.4.2 欧拉与龙格-库塔方法设有一阶方程与初始条件0(,)()y f x y y x y '=⎧⎨=⎩ (4.1) 其中(,)f x y 适当光滑,关于y 满足Lipschitz 条件,即存在L 使1212(,)(,)f x y f x y L y y -≤-则(4.1)式的解存在且惟一.关于()y y x =的解析解一般难以求到或根本无解析解,因而,实际问题中,通常,采用差分的方法. 在一系列离散点12n x x x <<<<上寻求其数值近似解12,,,,n y y y .相邻两个节点间的间距1n n n h x x +=-称为步长,一般地取等步长h ,则0,1,2,n x x nh n =+=1、欧拉方法在区间1[,]n n x x +上用差商1()()n n y x y x h+-代替(4.1)式中y ',对(,)f x y 中x 在1[,]n n x x +上取值n x 还是1n x +,而形成向前欧拉公式与向后欧拉公式.(1)向前欧拉公式(,)f x y 取左端点n x ,得如下公式1()()(,())n n n n y x y x hf x y x +≈+ (4.2)从0x 点出发,由初值00()y x y =代入(4.2)求得1000(,)y y hf x y =+ (4.3)反复利用(4.2),有1(,) 0,1,2,n n n n y y hf x y n +=+= (4.4)其几何意义如图4.1所示.图中()y y x =为方程(4.1)的精确 解曲线,其上任意点(,)x y 处切线斜率为(,)f x y . 从初值000(,)P x y 点出发,用该点斜率00(,)f x y 作一直线段,在1x x = 处得到111(,)P x y ,1y 由(4.2)式确定, 再从111(,)P x y 出发,以11(,)f x y 为斜率 作直线段,在2x x =处得到222(,)P x y ,y0yO0x 误差4x1x 2x 3x x()y y x =3()y x32()y y x - 0P1P2P3P 4P3y这一过程继续下去,形成折线012P PP ,作为积分曲线()y y x =的近似,用1()n y x + 图4.1 表示在1n x +处的精确值,1n y +为解的近似值,不难得到23211()()()()2n n n h y x y y x O h O h ++''-=+=这一误差称为局部截断误差. 若一种算法局部截断误差为1()P O h +,则称该算法具有P 阶精度,所以向前欧拉公式具有1阶精度.(2)向后欧拉公式若(,)f x y 中x 取1[,]n n x x +中的1n x +,则有如下公式:111(,) 0,1,2,n n n n y y hf x y n +++=+= (4.5)称式(4.5)为向后欧拉公式,因为此式中1n y +未知,故称其为隐式公式,无法用其直接计算1n y +,一般用向前欧拉公式产生初值.(0)11(,) 0,1,2,n n n n y y hf x y n ++=+=再按下式迭代(1)()111(,),0,1,,0,1,k k n n n n y y hf x y k n ++++=+==其误差估计如下23211()()()()2n n n h y x y y x o h o h ++''-=+=精度亦为1阶,将向前欧拉公式(4.4)与向后欧拉公式(4.5)及它们的误差的几何说明作一对比,是十分有益的,见图4.2. 为讨论局部截断误差,在图4.2中设点 (,)n n n P x y 落在积分曲线()y y x =上,按式 (4.4)及式(4.5)分别得1n P +点为A 与B , 且,A B 点一定在积分曲线()y y x =上相应 点θ的上、下两边,所以将式(4.4)与(4.5) 平均之,一定能得到更好的结果. (3)梯形公式 将向前与向后欧拉公式加以平均得到所 图4.2谓梯形公式111[(,)(,)] 0,1,2,2n n n n n n hy y f x y f x y n +++=++= (4.6)其局部截断误差为3()O h ,具有2阶精度.(4)改进的欧拉公式为使计算简单,又免去迭代的繁复,将公式(4.6)简化为两步1(,)n n n n y y hf x y +=+x 1n x +n x y()y y x =n PAθB111[(,)(,)], 0,1,2,2n n n n n n hy y f x y f x y n +++=++= (4.7)或写为1121211()2(,)(,)n n n n n n h y y k k k f x y k f x y hk ++⎧=++⎪⎪=⎨⎪=+⎪⎩0,1,2,n= (4.8)最后指出,上述欧拉方法可推广至微分方程组,如0000(,,)(,,)()()y f x y z z g x y z y x y z x z '=⎧⎪'=⎪⎨=⎪⎪=⎩ 向前欧拉公式为11(,,)(,,)n n n n n n n n n n y y hf x y z z z hg x y z ++=+⎧⎨=+⎩ 0,1,2,n=2、龙格_库塔方法 由微分中值定理1[()()]/(),01n n n y x y x h y x h θθ+'-=+<<又因为(,)y f x y '=,所以()(,())n n n y x h f x h y x h θθθ'+=++从而有1()()(),()n n n n y x y x hf x h y x h θθ+=+++ (4.9)令(,())n n k f x h y x h θθ=++,称其为区间1[,]n n x x +上的平均斜率,由(4.9)可知,给出一种平均斜率,可相应导出一种算法. 向前欧拉公式中(,)n n k f x y =,精度低. 改进欧拉公式中取111((,)(,))2n n n n k f x y f x y ++=+,精度提高,下面,我们在区间1[,]n n x x +内多取几个点,将其斜率加权平均,就能构造出精度更高的计算公式,公式的推导不再具体给出,只开列具体结果. (1)2阶龙格_库塔公式11122121()(,)(,),0,1n n n n n n y y h k k k f x y k f x ah y hk a λλββ+=++⎧⎪=⎨⎪=++<<⎩ (4.10) 其中12211,,12a aβλλλ+===,由于4个未知数只有3个方程,所以解不惟一,若令121,12a λλβ+===,即得改进的欧拉公式,具有2阶精度.(3)4阶龙格_库塔公式 只给出精典格式中最常用的一种.112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y k k k k k f x y h h k f x y k h h k f x y k k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ (4.11) 其计算精度为4阶4.3 模型与实验1、模型与问题 [例2] 单摆运动图4.3中一根长l 的细线,一端固定,另一端悬挂质量为m 的小球,在重力作用下,小球处于竖直的平衡位置. 现使小球偏离平衡位置一个小的角度θ,然后使其自由运动,在不考虑空气阻力情形下,小球将沿弧线作周期一定的简谐运动. 0θ=为平衡位置,在小球摆动过程中,当与平衡位置夹 角为θ时,小球所受重力在其动运轨迹的分量为sin mg θ- (负号表示力的方向使θ减少),由牛顿第二定律可得微分方 程()sin ml t mg θθ''=- (4.12) 设小球初始偏离角度为0θ,且初速为0,式(4.12)的初 始条件为0(0),(0)0θθθ'== (4.13)当0θ不大时,sin θθ≈,式(4.12)化为线性常系数微分方程 图4.30g lθθ''+= (4.14)解得θlθmg0()t θθ= (4.15)简谐运动的周期为2T =. 现在的问题是:当0θ较大时,仍用θ近似sin θ,误差太大,式(4.12)又无解析解,试用数值方法在0030,10θθ==两种情况下求解,画出()t θ的图形,与近似解(4.15)比较,这里设25cm l =. [例3] 捕食与被捕食当鲨鱼捕食小鱼,简记为乙捕食甲,在时刻t ,小鱼的数量为()x t ,鲨鱼的数量为()y t ,当甲独立生存时它的(相对)增长率与种群数量成正比,即有()()x t rx t '=,r 为增长率,而乙的存在使甲的增长率r 减少,设减少率与乙的数量成正比,而得微分方程()()(())x t x t r ay t rx axy '=-=- (4.16)比例系数a 反映捕食者掠取食饵的能力.乙离开甲无法生存,设乙独自存在时死亡率为d ,()()y t dy t '=-,甲为乙提供食物,使乙的死亡率d 降低,而促其数量增长,这一作用与甲的数量成正比,于是()y t 满足()(())y t y d bx t dy bxy '=-+=-+ (4.17)比例系数b 反映甲对乙的供养能力,设若甲,乙的初始数量分别为00(0);(0)x x y y == (4.18) 则微分方程(4.16),(4.17)及初始条件(4.18)确定了甲,乙数量()x t 、()y t 随时间变化而演变的过程,但该方程无解析解,试用数值解讨论以下问题:(1)设001,0.5,0.1,0.02,25,2r d a b x y ======,求方程(4.16),(4.17)在条件(4.18)下的数值解,画出(),()x t y t 的图形及相图(,)x y ,观察解(),()x t y t 的周期变化,近似确定解的周期和,x y 的最大、小值,近似计算,x y 在一个周期内的平均值. (2)从式(4.16)和(4.17)消去dt 得到()()dx x r ay dy y d bx -=-+ (4.19) 解方程(4.19),得到的解即为相轨线,说明这是封闭曲线,即解确为周期函数. (3)将方程(4.17)改写为1()[]y x t d b y'=+ (4.20)在一个周期内积分,得到()x t 一周期内的平均值,类似可得()y t 一周期内的平均值,将近似计算的结果与理论值比较. 2、实验(1)方程(单摆问题)0sin (0),(0)0ml mg θθθθθ''=-⎧⎨'==⎩无解析解,为求其数值解,先将其化为等价的一阶方程组. 令12,x x θθ'==,方程化为1221102sin (0),(0)0x x g x x l x x θ'=⎧⎪⎪'=-⎨⎪==⎪⎩ 其中,09.8,25,g l θ==为100.1745≈(弧度)与300.5236≈(弧度)两种情况,具体编程如下:先以danbai.m 为文件名存放待解方程. 键入: function xdot =danbai (t,x) %x=[x (1),x (2)],g=9.8;1=25; %x (1)为解向量, x (2)是导数 xdot (1)=x (2); %xdot (1)=x '(1)=θ' xdot (2)=-g/1*sin (x (1)); %xdot (2)=θ''xdot=[xdot (1);xdot (2)]; %必须将导数向量以列向量形式给出.再以主程序(求数值解)调用待求方程,主程序用main2.m 为文件名存盘,其代码如下: clear functions[t,x]=ode23 (‘danbai’,[0,10],[0.1745,0])%只计算0100.1745θ==(弧度)的情形. %对近似解0()cos ,t wt w θθ==2T = w=sqrt (9.8/25); y=0.1745*cos (w*t);[t,x (:,1),y] %显示数据,无分号. hold on %欲在同一幅图中画近似解. plot (t,x (:,1), ‘b’) %画数值解,绿色 plot (t,y, ‘g*’) %用*号,红色画近似解. Hold off(2)食饵_捕食者 方程()x t rx axy '=- ()y t dy bxy '=-+可化为如下形式00x r ay x d bx y y ⎡⎤-⎡⎤⎡⎤⎢⎥=⎢⎥⎢⎥⎢⎥-+⎣⎦⎣⎦⎣⎦初始条件00(0),(0)x x y y ==表示为00(0)(0)x x y y ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦以shier.m 存盘如下代码 function xdot=shier (t,x) r=1;d=0.5;a=0.1;b=0.02;xdot=diag ([r-a*x (2),-d+b*x (1)])*x;%x=(1)(2)x x x y ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦,xdot=x y ⎡⎤⎢⎥⎢⎥⎣⎦以main3.m 存盘如下代码. clear functions ts=0:0.1:15; x0=[25,2];[t,x]=ode45 (‘shier’,ts,x0);[t,x] %显示数据t,x,y plot (t,x)gird %加网格线 gtext (‘x(t)’),gtext (‘y(t)’), %用点鼠标方式pause,figure (2) %将x 1(t),x 2(t)放至指定点 plot (x (:,1),x (:,2)) %以x (1)与x (2)为坐标点画相图xlabel (‘x’),ylabel (‘y’)可以猜测,1()x t 与2()x t 是周期函数,相图12(,)x x 是封闭曲线,从数值解可近似定出周期为110.7,x 的最大和最小值分别为99.3与22.0,x 的最大和最小值分别为28.4和2.0,为求1x 与2x 在一个周期的平均值,只需键入: y1=x (1:108,1); %x 1周期为10.7. x1p=trapz (y1)*0.1/10.7, %trapz (y1)返回按 y2=x (1:108,2); %梯形法对y1的积 x2p=trapz (y2)*0.1/10.7, %分值10711()1(1)2i i y i y i x =++∆∑ 可得1225,10x x == 对方程()()dx x r ay dy y d bx -=-+化为d bx r aydx dy x y-+-=积分得解()()d bx r ay x e y e c --= (4.21)即为原方程组的相轨线,其中c 由初始条件确定. 为说明上述相轨线是封闭的,令();()d bx r ay f x x e g y y e --==设其最大值分别为,m m f g ,若00,x y 满足00(),()m m f x f f y g ==则有00,d rx y b a==(令0,0f g ''==,可解出00,x y ),又当m m f g c ⋅≥时,相轨线(4.21)有意义. 当m m f g c =时,相轨线退化为一个点00(,)P x y ,对0m m c f g <<时,相轨线如图 4.4(c ),而图(a ),(b )分别为()f x 与()g y 的图像,下面讨论如何由(a ),(b )画(c ).(a ) (b ) (c )图4.4设(0)m m c pg p f =<<,若令0y y =,则有()f x p =,由(a )知,12,x x ∃,使12()()f x f x p ==,且102x x x <<于是相轨线应过110(,)q x y 与220(,)q x y ,对12(,)x x x ∈,()f x p >,由()()()m m f x g x pg g y g =⇒<,令()g y q =,又由(b )知,存在12,y y 使12()()g y g y q ==,且102y y y <<,于是相轨线又过31(,)q x y 与42(,)q x y 两点,所以对12,q q 间每个x ,轨线总要通过纵坐标为12102,()y y y y y <<的两点,于是相轨线是一条封闭曲线.相轨线封闭等价于(),()x t y t 是周期函数,记周期为T ,为求其在一个周期内的平均值(,)x y ,由1()()yx t d b y=+两边在一个周期内积分有((0)())y y T =:011ln ()ln (0)()T y T y dT d x x t dt T T b b b-⎡⎤==+=⎢⎥⎣⎦⎰ 同理()f xmfP1x 0x 2xxyy 1y ()()n m n mf x fg y g ==2yyqm g()g vx1y 0y 2y1xx 2x1q 3q2q00(,)P x yr y a=从而 00,x x y y ==即(),()x t y t 的平均值恰为相轨线中心点p 的坐标,提请读者注意的是:这里的0x 与0y 与初始条件中的00,x y 不是一件事.注意到,,,r d a b 在生态学上的意义,上述结果表明,捕食者的数量与食饵的增长率成正比,与它掠取食饵的能力成反比,食饵的数量与捕食者死亡率成正比,与它供养捕食者的能力成反比.3、练习内容(1)编写改进欧拉公式求微分方程数值解的程序,并用其与ode23求下列微分方程数值解,对二者作出比较.a )22,(0)0y x y y '=-=或(0)1y =.b)222()0x y xy x n y '''++-=2()2,()22y y πππ'==-(Bessel 方程,这里令12n =,其精确解为sin y =. c)cos 0,(0)1,(0)0y y x y y '''+===. (2)倒圆锥形容器,上底面直径为1.2m ,容器的高亦为1.2m ,在锥尖的地方开有一直径为3cm 的小孔,容器装满水后,下方小孔开启,由水利学知识可知当水面高度为h 时,水从小孔中流出的速度为v g =为重力加速度,若孔口收缩系数为0.6(即若一个面积单位的小孔向外出水时,水柱截面积为0.6),问水从小孔中流完需多少时间?2分钟时,水面高度是多少?(3)一只小船渡过宽为d 的河流,目标是起点A 正对着的另一岸上B 点,已知河水流速1v 与船在静水中的速度2v 之比为k .(a )建立小船航线的方程,求其解析解.(b )设12100m,1m/s,2m/s d v v ===,用数值解法求渡河所需时间,任意时刻小船的位置及航行曲线,作图并与解析解比较.(c )若流速1v 为0,0.5,2 (m/s),结果将如何?(4)研究种群竞争模型. 当甲、乙两个种群各自生存时,数量演变服从下面规律1212()(1),()(1)x y x t r x y t r y n n =-=- 其中,(),()x t y t 分别为t 时刻甲,乙两个种群的数量,12,r r 为其固有增长率,12,n n 为它们的最大容量,而当这两个种群在同一环境中生存时,由于乙消耗有限资源而对甲的增长产生影响,将甲的方程修改为1112(1)x y x r x s n n =-- (4.22) 这里1s 的含意是:对于供养甲的资源而言,单位数量乙(相对2n )的消耗率为单位数量甲(相对1n )消耗的1s 倍,类似地,甲的存在亦影响乙的增长,乙的方程应改为2212()(1)x y y t r y s n n =-- (4.23) 给定种群的初始值为 00(0),(0)x x y y == (4.24)及参数121212,,,,,r r s s n n 后,方程(4.22)与(4.23)确定了两种群的变化规律,因其解析解不存在,试用数值解法研究以下问题:(a )设121212001,100,0.5,2,10r r n n s s x y ========,计算(),()x t y t ,画出它们的图形及相图(,)x y ,说明时间t 充分大以后(),()x t y t 的变化趋势(人们今天看到的已经是自然界长期演变的结果).(b )改变1200,,,r r x y ,但1s 与2s 不变(保持121,1s s <>),分析所得结果,若121.5(1),0.7(1)s s =>=<,再分析结果,由此你得到什么结论,请用各参数生态学上的含义作出解释.(c )试验当120.8(1),0.7(1)s s =<=<时会有什么结果;当121.5(1), 1.7(1)s s =>=>时又会出现什么结果,能解释这些结果吗?。

数值求解常微分方程PPT课件

数值求解常微分方程PPT课件

程叫做微分方程,求解微分方程必须附加某种
定解条件.微分方程和定解条件一起组成定解
问题,定解条件分为初始条件(初值问题)和边
界条件(边值问题)两种.未知函数为一元的微
分方程叫做常微分方程,未知函数为多元函数,
叫做偏微分方程.微分方程中导数的最高阶叫
做微分方程的阶.本章主要讨论一阶常微分方
程.
1
第1页/共51页
36
第36页/共51页
4. 阿达姆斯方法
我们已经知道,初值问题等价于积分方程, 即
y(xn1) y(xn )
xn1 xn
f (x, y(x))dx
对积分式分别采用矩形公式和梯形公式可得到 欧拉公式和改进的欧拉公式,截断误差分别为 O(h2)和O(h3)。
37
第37页/共51页
为此,我们自然可以想到,若用更高次的 插值多项式来代替f(x,y),则所得公式的精 度会更高。这就是线性多步法的起源思想。
本章前面介绍的方法称为单步法,因为 在计算yi+1时,只用到前面yi的值。而对于线 性多步法是要利用前面已经算出的若干个值yik,…,yi-1,yi来求yi+1。
38
第38页/共51页
现 用 k 次 多 项 式 Pk(x)
来代替f(x,y(x))
y(xi1) y(xi )
xi 1 xi
Pk (x)dx
为 了 改 善 精 度 , 将 函 数 y(x) 在 点 xi 处 的 导 数 y’(xi) 用 中 心 差 商 来 表示,即



式y变'
为( x:i误)

正y比( x于i h13), xi1

二y(阶x精i1度) xi1
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

f(xn,yn)
y y0 n 1 y(y x n 0 )h f(xn,yn) n0,1,2,
根据 y0 可以一步步计算出函数 y y(x) 在 x1, x2, x3 x4, …上的近似值 y1, y2, y3, y4 , …
常微分方程数值解是一组离散的函数值数据,它的 精确表达式很难求解得到,但可以进行插值计算后 用插值函数逼近 y(x)
四 常微分方程数值解法
1
整体概述
概述一
点击此处输入
相关文本内容
概述二
点击此处输入
相关文本内容
概述三
点击此处输入
相关文本内容
2
常微分方程数值解法
引言(常微分方程数值解法概述) 显式欧拉法、隐式欧拉法、二步欧拉法 局部截断误差与精度 改进的欧拉方法 龙格-库塔方法 收敛性与稳定性简述 一阶常微分方程组与高阶常微分方程
即积分区间为:[xn1, xn1],则:
xn1 xn1
ydxy(xn1)y(xn1)
xn1 xn1
f[x,y(x)]dx
(xn1xn1)f[xn,y(xn)] 中矩形公式
2hf[xn,y(xn)]
以 y(x) 在 xn 1, xn 上的近似值代替精确值可得:
yy0n1 y(yxn01)2hf(xn,yn)
3
引言
一阶常微分方程初值问题:
y f (x, y)
y
(
x0
)
y0
定理:若 f (x, y) 在某闭区域 R :
微分方程 初始条件
| x x 0 | a ,| y y 0 | b ( a 0 , b 0 )
上连续,且在 R 域内满足李普希兹 (Lipschitz) 条件, 即存在正数 L,使得对于 R 域内的任意பைடு நூலகம்值 y1, y2,下 列不等式成立:
➢ y(xn):待求函数 y(x) 在 xn 处的精确函数值 ➢ yn :待求函数 y(x) 在 xn 处的近似函数值
6
y ( x n ) l h i m 0 y ( x n h h ) y ( x n ) y ( x n 1 ) h y ( x n )
代入初值问题表达式可得:
yn1yn h
需要前两步 的 计 算 结11 果
梯形公式欧拉法
在数值积分法中,如果用梯形公式近似计算 f (x, y)
在区间 [xn, xn+1] 上的积分,即:
x x n n 1f[x,y(x)]d x(x n 1x n)f[x n,y(x n)] 2 f[x n 1,y(x n 1)]
h
2f[x n,y(x n)]f[x n 1,y(x n 1)]
梯形公式欧拉法: y n 1 y n h 2 [ ( y n x n 1 ) ( y n 1 x n 1 1 ) ]
[ 1 ( h 2 ) ] y n 1 [ 1 ( h 2 ) ] y n ( h 2 ) ( x n x n 1 2 )
y n 1 y n h ( y n 1 x n 1 1 )
( 1 h ) y n 1 y n h ( x n 1 1 ) y n h ( x n h 1 ) y n 1 ( 1 1 h ) [ y n h ( x n h 1 ) ] ( y n 0 .1 x n 0 .1 1 )1 .1
|f ( x , y 1 ) f ( x , y 2 ) | L |y 2 y 1 | 则上述初值问题的连续可微的解 y(x) 存在并且唯一。4
引言(续)
实际生产与科研中,除少数简单情况能获得初值问题 的初等解(用初等函数表示的解)外,绝大多数情况 下是求不出初等解的。 有些初值问题即便有初等解,也往往由于形式过于复 杂而不便处理。 实用的方法是在计算机上进行数值求解:即不直接求 y(x) 的显式解,而是在解所存在的区间上,求得一系 列点 xn (n 0, 1, 2, …) 上解的近似值。
7
欧拉方法(续)
方法二 数值积分法
将微分方程 y f (x, y) 在区间 [xn, xn+1] 上积分:
xn1 xn
ydxy(xn1)y(xn)
xn1 xn
f[x,y(x)]dx
(xn1xn)f[xn,y(xn)]
hf[xn,y(xn)]
同样以近似值 yn 代替精确值 y(xn) 可得:
(x n 1x n )f[x n 1 ,y (x n 1 )] h f[x n 1 ,y (x n 1 )]
这样便得到了隐式欧拉法: yn1 ynhf(xn1,yn1) y0 y(x0)
含有未知 的函数值
隐式欧拉法没有显式欧拉法方便
10
二步欧拉法
在数值积分法推导中,积分区间宽度选为两步步长,
yy0n1 y(yxn0) hf(xn,yn)
8
欧拉方法的几何意义
y P0
P5 P1 P2 P3 P4
P6 y y(x)
0 x0 x1 x2 x3 x4 x5 x6
x
9
隐式欧拉法
在数值积分法推导中,积分的近似值取为积分区间宽 度与右端点处的函数值乘积,即:
x x n n 1y d xy (x n 1 )y (x n )x x n n 1f[x ,y (x )]d x
用近似值代替精确值可得梯形公式欧拉法: h
y n 1 y n 2 [f(x n ,y n ) f(x n 1 ,y n 1 )]
上式右端出现了未知项,可见梯形法是隐式欧拉法
的一种;实际上,梯形公式欧拉法是显式欧拉法与
隐式欧拉法的算术平均。
12

用显式欧拉法、隐式欧拉法、梯形法求解初值问题:
y y x 1
5
欧拉(Euler)方法
方法一 化导数为差商的方法 y ( x n ) l h i m 0 y ( x n h h ) y ( x n ) y ( x n 1 ) h y ( x n )
由于在逐步求解的过程中,y(xn) 的准确值无法求解 出来,因此用其近似值代替。 为避免混淆,以下学习简记:
y(0)
1
取 h 0.1,计算到 x 0.5,并与精确解进行比较
解:由已知条件可得:h 0.1,x0 0, y0 1, f (x, y) y x 1
显式欧拉法:
yn1 yn h(yn xn 1)
(1h)yn hxn h
0.9yn 0.1xn 0.1
13
例:(续)
隐式欧拉法: 化简得:
相关文档
最新文档