如何求解时滞（时延）方程？

xsj3917

有这样子一个时滞方程（附图）：[dy(1)/dt;dy(2)/dt]=-C*[y(1);y(2)]+A*[tanh(y(1));tanh(y(2))]+B*[tanh(y(1)(t-T));tanh(y(2)(t-T)]
其中时滞T＝1+0.1*sin(t);
C=[1 0;0 1];
A=[2.0 -0.1;-5.0 3.0];
B=[-1.5 -0.1;-0.2 -2.5];
请教怎么样用MATLAB来求解呢？
是龙格－库塔法ode45吗？怎么写代码？
谢谢

[ 本帖最后由 xsj3917 于 2006-10-19 09:20 编辑 ]

fanghuikeer

ode45不能解你的方程，你参考一下matlab的帮助文件，你的方程要用dde23 去解
Solve delay differential equations (DDEs) with constant delays
我也没仔细看
自己加把劲，你的问题很简单的．．．

xsj3917

谢谢您的回复，dde23对应的时滞是常数，但是这里时滞项是一个变量，而非常量，所以说难办啊
能给个Example吗？

[ 本帖最后由 xsj3917 于 2006-10-19 10:37 编辑 ]

happy

参考matlab自带的例子ddex1

xsj3917

很遗憾ddex1 的时滞也是常数

happy

原帖由 xsj3917 于 2006-10-19 16:50 发表
很遗憾ddex1 的时滞也是常数

常数和非常数类似的，dde23的完整调用格式是
sol = dde23(ddefun,lags,history,tspan,options,p1,p2,...)

如果非常话，你可以通过p1,p2...进行传递

supermtd

对于时变时滞的微分方程需要用SIMLINK来做了，里面有专门的可变时滞的模块，楼主可以尝试一下。当然，不排除用M语言来实现，要麻烦些。但是DDE23还不知道如何解决时变时滞的问题。请6楼的朋友指点。

fandalei

龙格－库塔法不可以球吗？

haiyuntju

你好,我的课题中涉及到含有延时的微分代数方程的求解,我想请教一下在.m文件中如何求解这类方程呢?

mario_mm

你好，这个问题已经解决了吗?我现在也遇到了这个问题，麻烦看到留言回复，或者加我：83493430
谢谢！

xsj3917

因为Matlab运行速度慢的原因，我最终用C做的仿真，就是根据龙格－库塔法的公式用c语言写的代码

meiyongyuandeze

回复 11 # xsj3917 的帖子

可否分享下你的成果呢?

xsj3917

double *lgkt(double x, double x1, double y, double y1)  //龙格－库塔法迭代
{
        double* u=new double[2];
        kx1=-x+A[0][0]*g(x)+A[0][1]*g(y)+B[0][0]*g(x1)+B[0][1]*g(y1);
        ky1=-y+A[1][0]*g(x)+A[1][1]*g(y)+B[1][0]*g(x1)+B[1][1]*g(y1);
        kx2=-(x+h*kx1/2)+A[0][0]*g(x+h*kx1/2)+A[0][1]*g(y+h*ky1/2)
                +B[0][0]*g(x1+h*kx1/2)+B[0][1]*g(y1+h*ky1/2);
        ky2=-(y+h*ky1/2)+A[1][0]*g(x+h*kx1/2)+A[1][1]*g(y+h*ky1/2)
                +B[1][0]*g(x1+h*kx1/2)+B[1][1]*g(y1+h*ky1/2);
        kx3=-(x+h*kx2/2)+A[0][0]*g(x+h*kx2/2)+A[0][1]*g(y+h*ky2/2)
                +B[0][0]*g(x1+h*kx2/2)+B[0][1]*g(y1+h*ky2/2);
        ky3=-(y+h*ky2/2)+A[1][0]*g(x+h*kx2/2)+A[1][1]*g(y+h*ky2/2)
                +B[1][0]*g(x1+h*kx2/2)+B[1][1]*g(y1+h*ky2/2);
        kx4=-(x+h*kx3)+A[0][0]*g(x+h*kx3)+A[0][1]*g(y+h*ky3)
                +B[0][0]*g(x1+h*kx3)+B[0][1]*g(y1+h*ky3);
        ky4=-(y+h*ky3)+A[1][0]*g(x+h*kx3)+A[1][1]*g(y+h*ky3)
                +B[1][0]*g(x1+h*kx3)+B[1][1]*g(y1+h*ky3);
        x=x+h*(kx1+2*kx2+2*kx3+kx4)/6;
        y=y+h*(ky1+2*ky2+2*ky3+ky4)/6;
        u[0]=x;
        u[1]=y;
        return u;
}
g()是双曲正切函数

xunyuqi

暂时看不太懂，但是给了我方向