解微分方程组Word文档格式.docx
- 文档编号:15005828
- 上传时间:2022-10-26
- 格式:DOCX
- 页数:13
- 大小:233.59KB
解微分方程组Word文档格式.docx
《解微分方程组Word文档格式.docx》由会员分享,可在线阅读,更多相关《解微分方程组Word文档格式.docx(13页珍藏版)》请在冰豆网上搜索。
1.ode45求解的上手例子:
求解方程组
Dx=y+x(1-x^2-y^2);
Dy=-x+y*(1-x^2-y^2)
初值x=0.1;
y=0.2;
先说明一下最常用的ode45调用方式,和相应的函数文件定义格式。
[t,x]=ode45(odefun,tspan,x0);
其中,Fun就是导函数,tspan为求解的时间区间(或时间序列,如果采用时间序列,则必须单调),x0为初值。
这时,函数文件可以采用如下方式定义
functiondx=odefun(t,x)
对于上面的小例子,可以用如下的程序求解。
functionjixianhuan
clear;
clc
x0=[0.1;
0.2];
[t,x]=ode45(@jxhdot,[0,100],x0);
plot(x(:
1),x(:
2))
functiondx=jxhdot(t,x)
dx=[
x
(2)+x
(1).*(1-x
(1).^2-x
(2).^2);
-x
(1)+x
(2).*(1-x
(1).^2-x
(2).^2)
];
2.终值问题
tspan可以是递增序列,也可以为递减序列,若为递减则可求解终值问题。
[t,x]=ode45(@zhongzhiode,[3,0],[1;
0;
2]);
plot(t,x)
functiondx=zhongzhiode(t,x)
dx=[2*x
(2)^2-2;
-x
(1)+2*x
(2)*x(3)-1;
-2*x
(2)+2*x(3)^2-4];
结果如下
3.odeset
options=odeset('
name1'
value1,'
name2'
value2,...)
[t,x]=solver(@fun,tspan,x0,options)
通过odeset设置options
第一,通过求解选项的设置可以改善求解精度,使得原本可能不收敛的问题收敛。
options=odeset('
RelTol'
1e-10);
第二,求解形如M(t,x)x'
=f(t,x)的方程。
例如,方程
=-0.2x+yz+0.3xy
4.带附加参数的ode45
有时我们需要研究微分方程组中的参数对于解的影响,这时采用带有参数的ode45求解会使求解、配合循环使用,可以使得求解的过程更加简捷。
使用方法:
只需将附加参数放在options的后面就可以传递给odefun了。
看下面的例子。
functionRossler
a=[0.2,0.2];
b=[0.2,0.5];
c=[5.7,10];
x0=[000];
forjj=1:
2
[t,x]=ode45(@myRossler,[0,100],x0,[],a(jj),b(jj),c(jj));
figure;
plot3(x(:
2),x(:
3));
gridon;
end
functiondx=myRossler(t,x,a,b,c)
-x
(2)-x(3);
x
(1)+a*x
(2);
b+(x
(1)-c)*x(3)];
5.刚性方程的求解
刚性方程就是指各个自变量的变化率差异很大,会造成通常的求解方法失效。
这是matlab中自带的一个例子,使用ode15s求解,如果用ode45求解就会出现错误。
functionmyode15study
[t,Y]=ode15s(@vdp1000,[03000],[20]);
plot(T,Y(:
1),'
-o'
)
plot(Y(:
1),Y(:
functiondy=vdp1000(t,y)
dy=zeros(2,1);
dy
(1)=y
(2);
dy
(2)=1000*(1-y
(1)^2)*y
(2)-y
(1);
6.高阶微分方程的求解
通常的方法是进行变量替换,将原方程降阶,转换成更多变量的一阶方程组进行求解。
在这个例子里我们求解一个动力学系统里最常见的一个运动方程
其中f=sin(t)
functionmyhighoder
x0=zeros(6,1);
[t,x]=ode45(@myhigh,[0,100],x0);
plot(t,x(:
1))
functiondx=myhigh(t,x)
f=[sin(t);
0];
;
M=eye(3);
C=eye(3)*0.1;
K=eye(3)-0.5*diag(ones(2,1),1)-0.5*diag(ones(2,1),-1);
dx=[x(4:
6);
inv(M)*(f-C*x(4:
6)-K*x(1:
3))];
7.延迟微分方程
matlab提供了dde23求解非中性微分方程。
dde23的调用格式如下:
sol=dde23(ddefun,lags,history,tspan)
lags是延迟量,比如方程中包含y1(t-0.2)和y2(t-0.3)则可以使用lags=[0.2,0.3]。
这里的ddefun必须采用如下的定义方式:
dydt=ddefun(t,y,Z)
其中的Z(:
1)就是y(t-lags
(1)),Z(:
2)就是y(t-lags
(2))...
下面是个使用dde23求解延迟微分方程的例子。
functionmydde23study
%
Thedifferentialequations
%
y'
_1(t)=y_1(t-1)
_2(t)=y_1(t-1)+y_2(t-0.2)
_3(t)=y_2(t)
aresolvedon[0,5]withhistoryy_1(t)=1,y_2(t)=1,y_3(t)=1for
t<
=0.
lags=[1,0.2];
history=[1;
1;
1];
tspan=[0,5];
sol=dde23(@myddefun,lags,history,tspan)
plot(sol.x,sol.y)
functiondy=myddefun(t,y,Z)
dy=[
Z(1,1);
Z
(1)+Z(2,2);
y
(2)
8.ode15i求解隐式微分方程
[T,Y]=ode15i(odefun,tspan,y0,yp0)
yp0为y'
的初值。
odefun的格式如下
dy=odefun(t,y,yp),yp表示y'
,而方程中应该使得f(t,y,y'
)=0
functionmyodeIMP
Theproblemis
y
(1)'
=-0.04*y
(1)+1e4*y
(2)*y(3)
y
(2)'
=
0.04*y
(1)-1e4*y
(2)*y(3)-3e7*y
(2)^2
y(3)'
3e7*y
(2)^2
Itistobesolvedwithinitialconditionsy
(1)=1,y
(2)=0,y(3)=0
tosteadystate.
y0=[1;
fixed_y0=[1;
yp0=[000];
fixed_yp0=[];
[y0mod,yp0mod]=decic(@myodefunimp,0,y0,fixed_y0,yp0,fixed_yp0);
tspan=[0,logspace(-6,6)];
[t,y]=ode15i(@myodefunimp,tspan,y0mod,yp0mod);
y(:
2)=1e4*y(:
2);
semilogx(t,y)
functionres=myodefunimp(t,y,yp)
res=[
-yp
(1)-0.04*y
(1)+1e4*y
(2)*y(3);
-yp
(2)+0.04*y
(1)-1e4*y
(2)*y(3)-3e7*y
(2)^2;
-yp(3)+3e7*y
(2)^2;
这次要接触一个新的求解ode的方法,就是使用simulink的积分器求解。
1.还是做我们研究过的一个例子(在初识matlab微分方程
(2)中采用的)。
积分器中设置初始条件;
f(u)中指定Dx,Dy的计算公式。
运行这个仿真,scope中可以看到两个变量的时程如下:
在WorkSpace里可以得到tout和yout,执行plot(yout(:
1),yout(:
1))得到与ode45求解相似的结果如下
2.这部分解决一个使用ode求解器dde23没法求解的一类延迟微分方程(中性微分方程)。
形如x'
(t)=f(x'
(t-t1),x(t),x(t-t2),x(t-t3))这类方程。
dde23是无法求解的,但是可以借助simulink仿真求解。
看下面的这个例子。
(t)=A1*x(t-t1)+A2*x'
(t-t2)+B*u(t)
t1=0.15;
t2=0.5
A1=[-12
3
-3]
A2=[0.02
0
0]
B=[0]
[106
-116
62]
[0
0.03
[1]
[207
-207
113]
0.04]
[2]
在continuous里找到transportDelay,就可以实现对于信号的延迟,因此可以建立如下仿真模型
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分 方程组