使用 odeToVectorField 数值求解一对耦合的二阶 ODES

Posted 2023-02-18

【中文标题】使用 odeToVectorField 数值求解一对耦合的二阶 ODES

【英文标题】：Numerically solving a pair of coupled second order ODES with odeToVectorField

【发布时间】：2022-01-17 07:43:51

【问题描述】：

我正在尝试使用 MATLAB 中的一些函数来数值求解以下形式的一对耦合二阶 ODE

\ddotx = f(x,y,\dotx,\doty)

\ddoty = f(x,y,\dotx,\doty)。

我只能使用一个二阶 ODE，但我尝试编写的代码不适用于一对 ODE。

函数 odeToVectorField 有效地采用二阶 ODE，并将其写为一对耦合的一阶 ODE 的向量。 ode45 是通常的龙格-库塔解法。 xInit 和 yInit 对应于 x 和 y 的初始条件，然后目标是在特定时间间隔内绘制 x 和 y 对时间的曲线。

 gamma1=0.1;
    gamma2=0.1;
    a=1;
    m=1;
    g=9.8;
    d=1;

syms x(t) y(t)
eqn1=diff(x,2)== (gamma1*diff(x))/(a + m*d^2 + (m/2)*d^2*cos(y-x))  + (gamma2*diff(y))/(a+ (m/2)*cos(y-x))  -  ( (m/2)*d^2*sin(y-x)*(diff(x)^2 - diff(y)^2))/(a + m*d^2 + (m/2)*d^2*cos(y-x))    -    ((m/2)*d^2*diff(x)^2*(y-x))/(a+ (m/2)*cos(y-x))  -   ((m/2)*d*(3*g*sin(x) + g*sin(y)))/(a + m*d^2 + (m/2)*d^2*cos(y-x)) - ((m/2)*d*g*sin(y))/(a+ (m/2)*cos(y-x))              
     
eqn2=diff(y,2)==  (gamma1*diff(x))/((m/2)*d^2*cos(y-x)) +   (gamma2*diff(y))/a    -  ( (m/2)*d^2*sin(y-x)*(diff(x)^2 - diff(y)^2))/((m/2)*d^2*cos(y-x))   - ((m/2)*d^2*diff(x)^2*(y-x))/a    -  ((m/2)*d*(3*g*sin(x) + g*sin(y)))/((m/2)*d^2*cos(y-x)) -  ((m/2)*d*g*sin(y))/a           

V = odeToVectorField(eqn1,eqn2)

 M = matlabFunction(V,'vars','t','Y')           

  interval = [0 20];
  xInit = [2 0];
 yInit = [2 0];
   ySol = ode45(M,interval,xInit, yInit);            
tValues = linspace(0,20,100);
 yValues = deval(ySol,tValues,1);
 plot(tValues,yValues)

【问题讨论】：

当您可以轻松地为一阶系统编写函数时，为什么还要对odeToVectorField 使用符号操作？最后，这在数值求解器中可能会更快。

没有什么特别的原因，我只是觉得使用 MATLAB 已经有的内置函数可能更简洁、更容易？

【参考方案1】：

只是为了比较，不使用符号表达式，可以将这个等式实现为

function dV = M(t,V)
    gamma1=0.1;
    gamma2=0.1;
    a=1;
    m=1;
    g=9.8;
    d=1;
    
    x = V(1); dx = V(2); y = V(3); dy = V(4);
    ddx = (gamma1*dx)/(a + m*d^2 + (m/2)*d^2*cos(y-x))  + (gamma2*dy)/(a+ (m/2)*cos(y-x))  -  ( (m/2)*d^2*sin(y-x)*(dx^2 - dy^2))/(a + m*d^2 + (m/2)*d^2*cos(y-x))    -    ((m/2)*d^2*dx^2*(y-x))/(a+ (m/2)*cos(y-x))  -   ((m/2)*d*(3*g*sin(x) + g*sin(y)))/(a + m*d^2 + (m/2)*d^2*cos(y-x)) - ((m/2)*d*g*sin(y))/(a+ (m/2)*cos(y-x));          

    ddy =  (gamma1*dx)/((m/2)*d^2*cos(y-x)) +   (gamma2*dy)/a    -  ( (m/2)*d^2*sin(y-x)*(dx^2 - dy^2))/((m/2)*d^2*cos(y-x))   - ((m/2)*d^2*dx^2*(y-x))/a    -  ((m/2)*d*(3*g*sin(x) + g*sin(y)))/((m/2)*d^2*cos(y-x)) -  ((m/2)*d*g*sin(y))/a;
    dV = [dx ddx dy ddy];  
end%function

interval = [0 20];
xInit = [2 0];
yInit = [2 0];
vSol = ode45(M,interval,[ xInit yInit]);            
tValues = linspace(0,20,100);
xValues = deval(vSol,tValues,1);
plot(tValues,xValues)

这可行，但报告t=0.244 和t=0.245 之间存在奇点。