MatLa求解延迟微分方程的注意事项.docx
- 文档编号:9493535
- 上传时间:2023-05-19
- 格式:DOCX
- 页数:21
- 大小:52.60KB
MatLa求解延迟微分方程的注意事项.docx
《MatLa求解延迟微分方程的注意事项.docx》由会员分享,可在线阅读,更多相关《MatLa求解延迟微分方程的注意事项.docx(21页珍藏版)》请在冰点文库上搜索。
MatLa求解延迟微分方程的注意事项
MatLab求解延迟微分方程的注意事项
2010-08-1417:
04
MatLab求解延迟微分方程的注意事项
使用MATLAB求解延时微分方程的两种方法:
DDE23和SimuLink有些不同点需要注意,否则结果会出现错误
使用MATLAB来求解延迟微分方程是在生物数学和化学计算求解中经常遇到的事,在其它领域也比较常见。
我所知道的,在MATLAB中求解微分方程有三种方法:
1.使用ode45(龙格-库塔法的一个变种)求解,这时用一个数组,记录y的延迟项,但是c的值要考虑步长,再代入方程就能实现延时效应;2.使用dde23求解常数延时方程、使用ddesd可以用来求解延迟与时间t有关的延迟微分方程;3.使用SimuLink建模求解,SimuLink是求解广义微分代数系统的通用工具,功能很强大,但是看惯了编程指令的人可能不大习惯,调试似乎也不太方便。
既然本文专门讨论求解延迟微分方程,就先介绍一下专用工具dde23,dde系列求解函数是由SouthernMethodistUniversity的L.F.Shampine和S.Thompson根据他们早期使用Fortran编写的Fortran90DDESolver移植到MATLAB上的,从MATLAB6.5开始加入MATLAB的官方发行版,根据薛定宇教授在其几本关于MATLAB的著作中提到的,该函数返回的sol中结构体sol.x和sol.y均为按行排列,与ode45等不同,不太规范(没办法,因为这个函数本来就不是Mathworks的官方作品),不过这一点已经不大可能得到改进了,因为L.F.Shampine和S.Thompson已经决定停止维护这个文件。
如果您想进一步了解该函数,可以访问它的主页。
MATLAB帮助中关于该函数的介绍不很清楚,如果需要进一步了解这个函数,需要下载作者为其写的手册。
下面以MATLAB中所附的一个例程来说明这个函数与Simulink建模求解的不同。
在MATLABprompt中键入editddex1就会找看到函数作者所写的一个入门例子:
functionddex1
%DDEX1Example1forDDE23.
%ThisisasimpleexampleofWille'andBakerthatillustratesthe
%straightforwardformulation,computation,andplottingofthesolution
%ofasystemofdelaydifferentialequations(DDEs).
%
%Thedifferentialequations
%
%y'_1(t)=y_1(t-1)
%y'_2(t)=y_1(t-1)+y_2(t-0.2)
%y'_3(t)=y_2(t)
%
%aresolvedon[0,5]withhistoryy_1(t)=1,y_2(t)=1,y_3(t)=1for
%t<=0.
%
%Thelagsarespecifiedasavector[1,0.2],thedelaydifferential
%equationsarecodedinthesubfunctionDDEX1DE,andthehistoryis
%evaluatedbythefunctionDDEX1HIST.Becausethehistoryisconstantit
%couldbesuppliedasavector:
%sol=dde23(@ddex1de,[1,0.2],ones(3,1),[0,5]);
%
%SeealsoDDE23,FUNCTION_HANDLE.
%JacekKierzenka,LawrenceF.ShampineandSkipThompson
%Copyright1984-2004TheMathWorks,Inc.
%$Revision:
1.2.4.2$$Date:
2005/06/2119:
24:
16$
sol=dde23(@ddex1de,[1,0.2],@ddex1hist,[0,5]);
figure;
plot(sol.x,sol.y)
title('AnexampleofWille''andBaker.');
xlabel('timet');
ylabel('solutiony');
%--------------------------------------------------------------------------
functions=ddex1hist(t)
%ConstanthistoryfunctionforDDEX1.
s=ones(3,1);
%--------------------------------------------------------------------------
functiondydt=ddex1de(t,y,Z)
%DifferentialequationsfunctionforDDEX1.
ylag1=Z(:
1);
ylag2=Z(:
2);
dydt=[ylag1
(1)
ylag1
(1)+ylag2
(2)
y
(2)];
这里先不管函数使用的具体语法,求解模型为:
显然有两个延时常数1、0.2。
在使用手册中这样介绍:
theDDEreducestoaninitialvalueproblemforanODEwithy(t−1)equaltothegivenhistoryS(t−1)andinitialvaluey(0)=1
也就是说,仿真求解开始时间为t=0,此时y_1=1,要特别注意,这里设定的y的history并不仅仅是y的初始值,也包括t<0任意时刻的值。
也就是说,仿真开始时刻的
这样,仿真从0时刻开始y_1就以初始为1的变化率上升。
同理,y_2的情况也一样。
而在simulink中,我们搭建等效模型如图:
建模方程是正确的,各积分模块的初始值为1,延迟模块的延时值设为1和0.2,但是得出的结果与DDE23有非常大的差别。
经过数次更正求解器参数依然没有改观。
经过对比仿真结果可以发现:
不同于DDE23,SimuLink的仿真时间设为T=0,此时,y_1的初始值为1,但是Y_history的值没有得到体现,也就是说y'_1的值为0,只有到T=1时刻以后,才会有y'_1(t)=y(t-1),方程延时效应才正式得到体现。
这与DDE23的基本设定是有差别的。
如果方程模型中只有一个延时常数1,将simulink的仿真时间定为0-6秒,其中1-6秒与DDE23的仿真结果相同。
在本例中有两个时间常数,这样做行不通。
那么我们可以通过设定两个延时模块的initialoutput选项分别为y_1(-1)=1,y_2(-0.2)=1,这样就会得到与DDE23相同的结果。
因此,在使用MATLAB求解延迟微分方程时,应该先弄清楚延迟方程的具体意义,再选择适当的求解方式,才能得到正确的结果。
顺便提一下,MATLAB被广为诟病的除了昂贵的价格之外,最突出的就是透明度较低的内部算法,使用者往往搞不清算法的内部结构(很多算法甚至不予公开),也就无从判断求解结果的正确性,这时,转向一些开源软件(SciLab、MAXIMA、Octave、SAGE)和开放式求解库可以解决问题。
用dde求解延迟微分方程中
这是看到的资料:
functiondx=ddefun(t,x,z)
%z(i,j)表示状态变量xj延迟了tau(i)
tau1=z(:
1);%第一列表示延迟了tau
(1)=1的所有状态变量
tau2=z(:
2);
%x2延迟了1,故表示为z(2,1)或者tau1
(2)
%x1延迟了0.5,故表示为z(1,2)或者tau2
(1)
dx=[1-3*x
(1)-tau1
(2)-0.2*tau2
(1)^3-tau2
(1)
x(3)
4*x
(1)-2*x
(2)-3*x(3)];
z(i,j)表示状态变量xj延迟了tau(i),%x2延迟了1,故表示为z(2,1)。
这两句话不是矛盾吗?
因为tau=[10.5]
所以z(2,1)是不是应该表示x1延迟了0.5?
而z(1,2)表示的是x2延迟1?
而在同样求解同一个问题时:
functiondy=ddefun1(t,y,Z);
tau1=z(:
1);
d=0.1;b=0.002;a=5;p=0.05;c=0.2;b=0.3;
dy=[270-d*y
(1)-b*y
(1)*y
(2)
b*y
(1)*y
(2)-a*y
(2)-p*y
(2)*y(3)
c*tan1
(2)-b*y(3)];
end
functiondy=ddefun1(t,y,Z);
d=0.1;b=0.002;a=5;p=0.05;c=0.2;b=0.3;
dy=[270-d*y
(1)-b*y
(1)*y
(2)
b*y
(1)*y
(2)-a*y
(2)-p*y
(2)*y(3)
c*Z(2,1)-b*y(3)];
end
这两程序感觉上是一样的,但是画出来的图相差这么大
这是一个延迟微分方程;
MATLAB可以解这类延迟微分方程,但是是数值解法;所以需要之到一个初始条件
x(0)的值;
你能给出x(0)的值我可以帮你解
首先编写关于延迟函数的M文件;
functiondx=yanchi(t,x,z)
tau=z;%定义延迟时间
dx=x*(1-tau);%延迟函数
接下来命令求解
>>tau=1;%给定延迟时间
>>history=0;%初始值
>>tapan[0,10];%求解时间范围
>>sol=dde23(@yanchi,tau,history,tapan);%延迟问题求解
>>plot(sol.x,sol.y);%作图
下面附上了图片x(0)=0和x(0)=2的情况
显然初始值不同结果不同,这就是为什么需要初始值的情况
延迟微分方程是:
dx(t)/dt=0.2*x(t-17)/(1+x(t-17)^10)-0.1x(t)
初始条件是x(0)=1.2
当t<0时。
x(t)=0
要MATLAB编写的程序
多谢!
!
!
!
!
!
!
!
!
其他回答共1条
2009-4-2719:
36okhz|六级
%保存为jiang.m
functiondxdt=jiang(t,x)
dxdt=0.2*x.*(t-17)./(1+x.*(t-17).^10)-0.1*x;
%执行
[t,x]=ode45(@jiang,[-2020],1.2)
plot(t,x,'-o')
functionf=lorenz_ext(t,X)
%
%Lorenzequation
%
%dx/dt=SIGMA*(y-x)
%dy/dt=R*x-y-x*z
%dz/dt=x*y-BETA*z
%
%IndemorunSIGMA=10,R=28,BETA=8/3
%Initialconditions:
x(0)=0,y(0)=1,z(0)=0;
%Referencevaluesfort=10000:
%L_1=0.9022,L_2=0.0003,LE3=-14.5691
%
%See:
%K.Ramasubramanian,M.S.Sriram,"Acomparativestudyofcomputation
%ofLyapunovspectrawithdifferentalgorithms",PhysicaD139(2000)72-86.
%
%--------------------------------------------------------------------
%Copyright(C)2004,GovorukhinV.N.
%Valuesofparameters
a=10;
c=28;
b=8/3;
r=-0.3;
x=X
(1);y=X
(2);z=X(3);w=X(4);
Y=[X(5),X(9),X(13),X(17);
X(6),X(10),X(14),X(18);
X(7),X(11),X(15),X(19);
X(8),X(12),X(16),X(20)];
f=zeros(16,1);
%Lorenzequation
f
(1)=a*(y-x)+w;
f
(2)=-x*z+c*x-y;
f(3)=x*y-b*z;
f(4)=-y*z+r*w;
%Linearizedsystem
Jac=[-a,a,0,1;
-z+c,-1,-x,0;
y,x,-b,0;
0,-z,-y,r];
%Variationalequation
f(5:
20)=Jac*Y;
clearY;
clearJac;
%Outputdatamustbeacolumnvector
function[Texp,Lexp]=lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp);
%
%LyapunovexponentcalcullationforODE-system.
%
%Thealogrithmemployedinthism-filefordeterminingLyapunov
%exponentswasproposedin
%
%A.Wolf,J.B.Swift,H.L.Swinney,andJ.A.Vastano,
%"DeterminingLyapunovExponentsfromaTimeSeries,"PhysicaD,
%Vol.16,pp.285-317,1985.
%
%ForintegratingODEsystemcanbeusedanyMATLABODE-suitemethods.
%ThisfunctionisapartofMATDSprogram-toolboxfordynamicalsysteminvestigation
%See:
http:
//www.math.rsu.ru/mexmat/kvm/matds/
%
%Inputparameters:
%n-numberofequation
%rhs_ext_fcn-handleoffunctionwithrighthandsideofextendedODE-system.
%ThisfunctionmustincludeRHSofODE-systemcoupledwith
%variationalequation(nitemsoflinearizedsystems,seeExample).
%fcn_integrator-handleofODEintegratorfunction,forexample:
@ode45
%tstart-startvaluesofindependentvalue(timet)
%stept-stepont-variableforGram-Schmidtrenormalizationprocedure.
%tend-finishvalueoftime
%ystart-startpointoftrajectoryofODEsystem.
%ioutp-stepofprinttoMATLABmainwindow.ioutp==0-noprint,
%ifioutp>0theneachioutp-thpointwillbeprint.
%
%Outputparameters:
%Texp-timevalues
%Lexp-Lyapunovexponentstoeachtimevalue.
%
%UsershavetowritetheirownODEfunctionsfortheirspecified
%systemsandusehandleofthisfunctionasrhs_ext_fcn-parameter.
%
%Example.Lorenzsystem:
%dx/dt=sigma*(y-x)=f1
%dy/dt=r*x-y-x*z=f2
%dz/dt=x*y-b*z=f3
%
%TheJacobianofsystem:
%|-sigmasigma0|
%J=|r-z-1-x|
%|yx-b|
%
%Then,thevariationalequationhasaform:
%
%F=J*Y
%whereYisasquarematrixwiththesamedimensionasJ.
%Correspondingm-file:
%functionf=lorenz_ext(t,X)
%SIGMA=10;R=28;BETA=8/3;
%x=X
(1);y=X
(2);z=X(3);
%
%Y=[X(4),X(7),X(10);
%X(5),X(8),X(11);
%X(6),X(9),X(12)];
%f=zeros(9,1);
%f
(1)=SIGMA*(y-x);f
(2)=-x*z+R*x-y;f(3)=x*y-BETA*z;
%
%Jac=[-SIGMA,SIGMA,0;R-z,-1,-x;y,x,-BETA];
%
%f(4:
12)=Jac*Y;
%
%RunLyapunovexponentcalculation:
%
%[T,Res]=lyapunov(3,@lorenz_ext,@ode45,0,0.5,200,[010],10);
%
%Seefiles:
lorenz_ext,run_lyap.
%
%--------------------------------------------------------------------
%Copyright(C)2004,GovorukhinV.N.
%ThisfileisintendedforusewithMATLABandwasproducedforMATDS-program
%http:
//www.math.rsu.ru/mexmat/kvm/matds/
%lyapunov.misfreesoftware.lyapunov.misdistributedinthehopethatit
%willbeuseful,butWITHOUTANYWARRANTY.
%
%
%n=numberofnonlinearodes
%n2=n*(n+1)=totalnumberofodes
%
n1=n;n2=n1*(n1+1);
%Numberofsteps
nit=round((tend-tstart)/stept);
%Memoryallocation
y=zeros(n2,1);cum=zeros(n1,1);y0=y;
gsc=cum;znorm=cum;
%Initialvalues
y(1:
n)=ystart(:
);
fori=1:
n1y((n1+1)*i)=1.0;end;
t=tstart;
%Mainloop
[TT,YY]=feval(fcn_integrator,rhs_ext_fcn,[tt+5000],y);
y=YY(size(YY,1),:
);
fprintf('yvalue');
fori=1:
n1
fprintf('%10.6f',Y(i));
end;
fprintf('\n');
forITERLYAP=1:
nit
%SolutuionofextendedODEsystem
[T,Y]=feval(fcn_integrator,rhs_ext_fcn,[tt+stept],y);
%if(mod(ITERLYAP,ioutp)==0)
%fprintf('yvalue');
%fori=1:
n1
%fprintf('%10.6f',Y(i));
%end;
%fprintf('\n');
%end;
%forii=1:
500
%ddy=lorenz_ext(1,y);
%Y=ddy*0.001+y;
%y=Y;
%end;
t=t+stept;
y=Y(size(Y,1),:
);
%y=Y;
fori=1:
n1
forj=1:
n1y0(n1*i+j)=y(n1*j+i);end;
end;
%
%constructneworthonormalbasisbygram-schmidt
%
znorm
(1)=0.0;
forj=1:
n1znorm
(1)=znorm
(1)+y0(n1*j+1)^2;end;
znorm
(1)=sqrt(znorm
(1));
forj=1:
n1y0(n1*j+1)=y0(n1*j+1)/znorm
(1);end;
forj=2:
n1
fork
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MatLa 求解 延迟 微分方程 注意事项