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
升级会员 