I am attempting to use some of the functions in MATLAB to numerically solve a pair of coupled second order ODEs of the form

\ddot{x} = f(x,y,\dot{x},\dot{y})

\ddot{y} = f(x,y,\dot{x},\dot{y}).

I am able to get it to work with just one second-order ODE, but the code I am trying to does not work for a pair of ODEs.

The function odeToVectorField effectively takes a second order ODE and writes it as a vector for a pair of coupled first order ODEs. ode45 is the usual Runge-Kutta solution method. xInit and yInit correspond to the initial conditions for x and y and the aim is then to plot both x and y against time over a certain interval of time.

 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)

CodePudding user response：

Just to compare, without using symbolic expressions one would implement this equation as

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)

This works, but reports a singularity between t=0.244 and t=0.245.