常微分方程数值解
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
30 25
rx axy x dy bxy y x ( 0) x , y ( 0) y 0 0
100 80
shier.m, shier1.m
x(t)
60
y
20 15 10
40 20 0
y(t)
0 5 10 15
5 0 0 20 40 x 60 80 100
猜测 x(t), y(t)是周期函数; y(x)是封闭曲线
数值积分计算一个周期的平均值： x 25, y 10
15
function y=shier(t,x) r=1;d=0.5;a=0.1;b=0.02; %y=diag([r-a*x(2),-d+b*x(1)])*x; y=[(r-a*x(2))*x(1),(-d+b*x(1))*x(2)]';
14.9986
y( t ) 0 0.5058 2.1308 4.8283 8.2755 12.0830
40.0164
y1 ( t ) 0 3.5 7.0 10.5 14.0 17.5
42.0
t=? 缉私艇追上走私船
累积误差较大
提高精度!
1.3 1.4 1.5 1.6
14.9996 15.0117 15.0023 14.9866
16
ts=0:0.1:15; x0=[25,2]; [t,x]=ode45('shier',ts,x0); [t,x], plot(t,x),grid, gtext('x(t)','FontSize',16),gtext('y(t)','FontSize',16), pause, plot(x(:,1),x(:,2)),grid, xlabel('x','FontSize',16),ylabel('y','FontSize',16) % averge value of x: T is about 10.7 %trapz(t(1:108),x(1:108,:))/10.7
2
模型的解析解
1 c x k c x k p [( ) ( ) ] 2 c c
p(0) 0
k a /b 1
y(0) 0
c 1 c x 1 k 1 c x 1 k kc y [ ( ) ( ) ] 2 1 k c 1 k c 1 k 2
缉私艇的航线y(x)的解析解
令 p
dy dx
dy dy dt dx dt dx
dp a 2 (c x ) k 1 p , k dx b p(0) 0
11
dp (c x ) k 1 p2 dx
实例1
海上缉私(续)
c x k 1 p p ( ) c cx k 2 1 p p ( ) c
7
实例1
海上缉私(续) 模型的数值解
t 0 0.05 0.10 0.15 0.20 0.25 0.30
x
a=20, b=40, c=15
走私船的位置 x1(t)= c=15 y1(t)=at=20t 缉私艇的航线y(x)
10 8 6 4 2 0 y
x(t) 0 1.9984 3.9854 5.9445 7.8515 9.6705 11.3496
17
实例2 弱肉强食
(r ay) x x (d bx) y y
弱肉 种群甲靠丰富的自然资源生存 食饵(Prey) 强食 种群乙靠捕食种群甲为生 捕食者(Predator) 两个种群的数量如何演变?
• 历史背景——一次世界大战期间地中海渔业 的捕捞量下降(食用鱼和鲨鱼同时捕捞)，但 是其中鲨鱼的比例却增加，为什么？
13
实 例 2
模型
/ x r, r 0 x
b(c x ) (c x )2 ( at y ) 2 b( at y ) (c x ) ( at y )
2 2
y Q(c,at) R(c,y ) c x
3
P(x,y) 0
b
由方程无法得到x(t), y(t)的解析解 需要用数值解法求解
龙格—库塔方法的 MATLAB 实现
常微分方程数值解
1
实例1 海上缉私
海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船 正以速度a向正北方向行驶，缉私艇立即以最大速度b(>a)前往 拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终 指向走私船。 • 建立任意时刻缉私艇位置及 航线的数学模型,并求解; • 求出缉私艇追上走私船的时间。
( r ay ) x rx axy x ( d bx ) y dy bxy y x(0) x , y (0) y 0 0
Volterra模型
x(t), y(t)无解析解
14
实例2
弱肉强食 模型的数值解
r 1, d 0.5 a 0.1, b 0.02 x0 25, y0 2
(t ) f (t , x), x(t0 ) x0 , x ( x1 , xn )T , f ( f1 , f n )T x
[t,x]=ode23(@f,ts,x0,opt)
3级2阶龙格-库塔公式
[t,x]=ode45(@f,ts,x0,opt) 5级4阶龙格-库塔公式 function dx=f(t,x) f是待解方程写成 dx=[f1; f2;; fn]; 的函数m文件： ts = [t0,t1, …,tf] 输出指定时刻 t0,t1, …,tf 的函数值 ts = t0:k:tf 输出 [t0,tf] 内等分点处的函数值
北 b a
艇
c
船
2
实例1 海上缉私
建立坐标系如图: t=0 艇在(0, 0), 船在(c, 0); 船速a, 艇速b 时刻 t 艇位于P(x, y), 船到达 Q(c, at)
模型: dx
dx dt dy dt
dy b cos , b sin dt dt
kc abc x=c时 y 2 2 2 缉私艇追上走私船的y坐标 1 k b a
缉私艇追上走私船的时间:
a=20, b=40, c=15 t1=0.5
t1
bc b2 a 2
a=35, b=40, c=15 t1=1.6
12
实 例 2
弱 肉 强 食
问 自然界中同一环境下两个种群之间的生存方式 题 相互竞争 相互依存 弱肉强食
44.000005 48.000005 52.000005 55.999931
45.5 49.0 52.5 56.0
10
实例1 海上缉私(续)
dx dt b (c x ) (c x) (at y )
2 2
模型的解析解
,
dy dt
b(at y ) (c x) 2 (at y ) 2
44.0165 48.0183 52.0146 55.9486
45.5 49.0 52.5 56.0
9
模型的数值解 t x(t) y( t ) 0 0 0 opt=odeset('RelTol',1e-6, 0.1 3.956104 0.505813 'AbsTol',1e-9); 0.2 7.592822 2.130678 [t,x]=ode45(@jisi,ts,x0,opt); 0.3 10.521921 4.829308 a=35, b=40, c=15 0.4 12.539454 8.269840 缉私艇的航线y(x) 60 y 0.5 13.753974 12.075344
6
• • • • • • • • • • • • • • • • • • • • •
clear,clf,shg %Set the definied time %ts=0:0.05:0.5; %ts=0:0.1:1.6; n=length(ts); x0=[0 0]; a=35;b=40;c=15; opt1=odeset('RelTol',1e-6,'AbsTol',1e-9); [t,x]=ode45(@jisi,ts,x0,opt1,a,b,c); %a=35;b=40;c=15; %opt1=odeset('RelTol',1e-6,'AbsTol',1e-9); %[t,x]=ode45(@jisi,ts,x0,opt1,a,b,c); %exact solution x1=c y1=a*t; %output t,x(t),y(t) and draw x(t),y(t) [t,x,y1] plot(t,x),grid,gtext('x(t)','FontSize',16), gtext('y(t)','FontSize',16),pause %draw y(x): the position of tatch jisi plot(x(:,1),x(:,2),'r*'),grid xlabel('x','FontSize',16),ylabel('y','FontSize',16)
x0为函数初值(n维) 输出t=ts, x为相应函数值(n维)
opt为选项，缺省时精度为：相对误差10-3， 绝对误差10-6, 计算步长按精度要求自动调整
4
实例1 海上缉私(续)
模型的数值解
dx dt dy dt b(c x ) (c x )2 ( at y ) 2 b( at y ) (c x )2 ( at y ) 2
y( t ) 0 0.0698 0.2924 0.6906 1.2899 2.1178 3.2005
y1 ( t ) 0 1.0 2.0 3.0 4.0 5.0 6.0
0
5
10
15
20
t=0.5时缉私艇追上走私船
0.35 0.40 0.45 0.50
12.8170 13.9806 14.7451 15.0046
y(0) 0
0 (x(t),y(t)) a b c x
y
x(0) 0,
设: 船速a=20 (海里/小时) 艇速b=40 (海里/小时) 距离c=15 (海里)
求: 缉私艇的位置x(t),y(t) 缉私艇的航线y(x)
jisi.m, seajisi.m
5
%Creat the function for jisi %Let x(1)=x, x(2)=y function dx=jisi(t,x,a,b,c) s=sqrt((c-x(1))^2+(a*t-x(2))^2); %dx=[b*(c-x(1))/s;b*(a*t-x(2))/s]; dx=[1;1]; dx(1)=b*(c-x(1))/s; dx(2)=b*(a*t-x(2))/s;
d2y dt (c x ) 2 a dx dx
dy at y dx (c x )
ds b dt
ds ( dx )2 ( dy ) 2
dy (c x) y at dx
dt dt ds dx ds dy
d2y a dy (c x ) 2 1 ( )2 dx b dx
4.5552 6.1773 8.0273 9.9979
7.0 8.0 9.0 10.0
8
实例1 海上缉私(续)
设b，c不变，a变大 为30, 35, 接近40, 观察解的变化 : a=35, b=40, c=15
模型的数值解
t 0 0.1 0.2 0.3 0.4 0.5
1.2
x(t) 0 3.9561 7.5928 10.5240 12.5384 13.7551
/ x r ay, a 0 x
弱 肉 强 食
甲独立生存的增长率r 乙使甲的增长率减小, 减小量与 y 成正比 乙独立生存的死亡率d
食饵(甲) 的密度x(t), 捕食者(乙)的密度y(t)
/ y d , d 0 y
/ y (d bx), b 0 甲使乙的死亡率减小, y 减小量与 x成正比
40
实例1 海上缉私(续)
y1 ( t ) 0 3.5 7.0 10.5 14.0 17来自百度文库5
20 x 0 0 5 10 15 20
1.2 14.999616 40.000005 42.0
t=1.6时缉私艇追上走私船
判断“追上”的有效方 法?
1.3 1.4 1.5 1.6
14.999963 14.999993 14.999998 15.000020
rx axy x dy bxy y x ( 0) x , y ( 0) y 0 0
100 80
shier.m, shier1.m
x(t)
60
y
20 15 10
40 20 0
y(t)
0 5 10 15
5 0 0 20 40 x 60 80 100
猜测 x(t), y(t)是周期函数; y(x)是封闭曲线
数值积分计算一个周期的平均值： x 25, y 10
15
function y=shier(t,x) r=1;d=0.5;a=0.1;b=0.02; %y=diag([r-a*x(2),-d+b*x(1)])*x; y=[(r-a*x(2))*x(1),(-d+b*x(1))*x(2)]';
14.9986
y( t ) 0 0.5058 2.1308 4.8283 8.2755 12.0830
40.0164
y1 ( t ) 0 3.5 7.0 10.5 14.0 17.5
42.0
t=? 缉私艇追上走私船
累积误差较大
提高精度!
1.3 1.4 1.5 1.6
14.9996 15.0117 15.0023 14.9866
16
ts=0:0.1:15; x0=[25,2]; [t,x]=ode45('shier',ts,x0); [t,x], plot(t,x),grid, gtext('x(t)','FontSize',16),gtext('y(t)','FontSize',16), pause, plot(x(:,1),x(:,2)),grid, xlabel('x','FontSize',16),ylabel('y','FontSize',16) % averge value of x: T is about 10.7 %trapz(t(1:108),x(1:108,:))/10.7
2
模型的解析解
1 c x k c x k p [( ) ( ) ] 2 c c
p(0) 0
k a /b 1
y(0) 0
c 1 c x 1 k 1 c x 1 k kc y [ ( ) ( ) ] 2 1 k c 1 k c 1 k 2
缉私艇的航线y(x)的解析解
令 p
dy dx
dy dy dt dx dt dx
dp a 2 (c x ) k 1 p , k dx b p(0) 0
11
dp (c x ) k 1 p2 dx
实例1
海上缉私(续)
c x k 1 p p ( ) c cx k 2 1 p p ( ) c
7
实例1
海上缉私(续) 模型的数值解
t 0 0.05 0.10 0.15 0.20 0.25 0.30
x
a=20, b=40, c=15
走私船的位置 x1(t)= c=15 y1(t)=at=20t 缉私艇的航线y(x)
10 8 6 4 2 0 y
x(t) 0 1.9984 3.9854 5.9445 7.8515 9.6705 11.3496
17
实例2 弱肉强食
(r ay) x x (d bx) y y
弱肉 种群甲靠丰富的自然资源生存 食饵(Prey) 强食 种群乙靠捕食种群甲为生 捕食者(Predator) 两个种群的数量如何演变?
• 历史背景——一次世界大战期间地中海渔业 的捕捞量下降(食用鱼和鲨鱼同时捕捞)，但 是其中鲨鱼的比例却增加，为什么？
13
实 例 2
模型
/ x r, r 0 x
b(c x ) (c x )2 ( at y ) 2 b( at y ) (c x ) ( at y )
2 2
y Q(c,at) R(c,y ) c x
3
P(x,y) 0
b
由方程无法得到x(t), y(t)的解析解 需要用数值解法求解
龙格—库塔方法的 MATLAB 实现
常微分方程数值解
1
实例1 海上缉私
海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船 正以速度a向正北方向行驶，缉私艇立即以最大速度b(>a)前往 拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终 指向走私船。 • 建立任意时刻缉私艇位置及 航线的数学模型,并求解; • 求出缉私艇追上走私船的时间。
( r ay ) x rx axy x ( d bx ) y dy bxy y x(0) x , y (0) y 0 0
Volterra模型
x(t), y(t)无解析解
14
实例2
弱肉强食 模型的数值解
r 1, d 0.5 a 0.1, b 0.02 x0 25, y0 2
(t ) f (t , x), x(t0 ) x0 , x ( x1 , xn )T , f ( f1 , f n )T x
[t,x]=ode23(@f,ts,x0,opt)
3级2阶龙格-库塔公式
[t,x]=ode45(@f,ts,x0,opt) 5级4阶龙格-库塔公式 function dx=f(t,x) f是待解方程写成 dx=[f1; f2;; fn]; 的函数m文件： ts = [t0,t1, …,tf] 输出指定时刻 t0,t1, …,tf 的函数值 ts = t0:k:tf 输出 [t0,tf] 内等分点处的函数值
北 b a
艇
c
船
2
实例1 海上缉私
建立坐标系如图: t=0 艇在(0, 0), 船在(c, 0); 船速a, 艇速b 时刻 t 艇位于P(x, y), 船到达 Q(c, at)
模型: dx
dx dt dy dt
dy b cos , b sin dt dt
kc abc x=c时 y 2 2 2 缉私艇追上走私船的y坐标 1 k b a
缉私艇追上走私船的时间:
a=20, b=40, c=15 t1=0.5
t1
bc b2 a 2
a=35, b=40, c=15 t1=1.6
12
实 例 2
弱 肉 强 食
问 自然界中同一环境下两个种群之间的生存方式 题 相互竞争 相互依存 弱肉强食
44.000005 48.000005 52.000005 55.999931
45.5 49.0 52.5 56.0
10
实例1 海上缉私(续)
dx dt b (c x ) (c x) (at y )
2 2
模型的解析解
,
dy dt
b(at y ) (c x) 2 (at y ) 2
44.0165 48.0183 52.0146 55.9486
45.5 49.0 52.5 56.0
9
模型的数值解 t x(t) y( t ) 0 0 0 opt=odeset('RelTol',1e-6, 0.1 3.956104 0.505813 'AbsTol',1e-9); 0.2 7.592822 2.130678 [t,x]=ode45(@jisi,ts,x0,opt); 0.3 10.521921 4.829308 a=35, b=40, c=15 0.4 12.539454 8.269840 缉私艇的航线y(x) 60 y 0.5 13.753974 12.075344
6
• • • • • • • • • • • • • • • • • • • • •
clear,clf,shg %Set the definied time %ts=0:0.05:0.5; %ts=0:0.1:1.6; n=length(ts); x0=[0 0]; a=35;b=40;c=15; opt1=odeset('RelTol',1e-6,'AbsTol',1e-9); [t,x]=ode45(@jisi,ts,x0,opt1,a,b,c); %a=35;b=40;c=15; %opt1=odeset('RelTol',1e-6,'AbsTol',1e-9); %[t,x]=ode45(@jisi,ts,x0,opt1,a,b,c); %exact solution x1=c y1=a*t; %output t,x(t),y(t) and draw x(t),y(t) [t,x,y1] plot(t,x),grid,gtext('x(t)','FontSize',16), gtext('y(t)','FontSize',16),pause %draw y(x): the position of tatch jisi plot(x(:,1),x(:,2),'r*'),grid xlabel('x','FontSize',16),ylabel('y','FontSize',16)
x0为函数初值(n维) 输出t=ts, x为相应函数值(n维)
opt为选项，缺省时精度为：相对误差10-3， 绝对误差10-6, 计算步长按精度要求自动调整
4
实例1 海上缉私(续)
模型的数值解
dx dt dy dt b(c x ) (c x )2 ( at y ) 2 b( at y ) (c x )2 ( at y ) 2
y( t ) 0 0.0698 0.2924 0.6906 1.2899 2.1178 3.2005
y1 ( t ) 0 1.0 2.0 3.0 4.0 5.0 6.0
0
5
10
15
20
t=0.5时缉私艇追上走私船
0.35 0.40 0.45 0.50
12.8170 13.9806 14.7451 15.0046
y(0) 0
0 (x(t),y(t)) a b c x
y
x(0) 0,
设: 船速a=20 (海里/小时) 艇速b=40 (海里/小时) 距离c=15 (海里)
求: 缉私艇的位置x(t),y(t) 缉私艇的航线y(x)
jisi.m, seajisi.m
5
%Creat the function for jisi %Let x(1)=x, x(2)=y function dx=jisi(t,x,a,b,c) s=sqrt((c-x(1))^2+(a*t-x(2))^2); %dx=[b*(c-x(1))/s;b*(a*t-x(2))/s]; dx=[1;1]; dx(1)=b*(c-x(1))/s; dx(2)=b*(a*t-x(2))/s;
d2y dt (c x ) 2 a dx dx
dy at y dx (c x )
ds b dt
ds ( dx )2 ( dy ) 2
dy (c x) y at dx
dt dt ds dx ds dy
d2y a dy (c x ) 2 1 ( )2 dx b dx
4.5552 6.1773 8.0273 9.9979
7.0 8.0 9.0 10.0
8
实例1 海上缉私(续)
设b，c不变，a变大 为30, 35, 接近40, 观察解的变化 : a=35, b=40, c=15
模型的数值解
t 0 0.1 0.2 0.3 0.4 0.5
1.2
x(t) 0 3.9561 7.5928 10.5240 12.5384 13.7551
/ x r ay, a 0 x
弱 肉 强 食
甲独立生存的增长率r 乙使甲的增长率减小, 减小量与 y 成正比 乙独立生存的死亡率d
食饵(甲) 的密度x(t), 捕食者(乙)的密度y(t)
/ y d , d 0 y
/ y (d bx), b 0 甲使乙的死亡率减小, y 减小量与 x成正比
40
实例1 海上缉私(续)
y1 ( t ) 0 3.5 7.0 10.5 14.0 17来自百度文库5
20 x 0 0 5 10 15 20
1.2 14.999616 40.000005 42.0
t=1.6时缉私艇追上走私船
判断“追上”的有效方 法?
1.3 1.4 1.5 1.6
14.999963 14.999993 14.999998 15.000020