Numerically solving a pair of coupled second order ODES with odeToVectorField

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 and also x against y over 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)

Answer 1

[Star Strider]

‘... the code I am trying to does not work for a pair of ODEs.’

Yes, it does. Look at the ‘Subs’ result from odeToVectorField to see that everything is there.

The problem that remains is that this is now a degree system so ‘xInit, yInit’ need to be concatenated with square brackets for it to work —

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

sympref('AbbreviateOutput',false);

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))              

eqn1(t) = 

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           

eqn2(t) = 

[V,Subs] = odeToVectorField(eqn1,eqn2)

V = 

Subs = 

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

M = function_handle with value:
    @(t,Y)[Y(2);sin(Y(1)).*(-4.9e+1./1.0e+1)+Y(4)./(cos(Y(1)-Y(3)).*5.0)-((Y(1)-Y(3)).*Y(4).^2)./2.0+Y(2)./1.0e+1-((sin(Y(1)).*(4.9e+1./1.0e+1)+sin(Y(3)).*(1.47e+2./1.0e+1)).*2.0)./cos(Y(1)-Y(3))+(sin(Y(1)-Y(3)).*(Y(2).^2-Y(4).^2))./cos(Y(1)-Y(3));Y(4);(sin(Y(1)).*(-4.9e+1./1.0e+1))./(cos(Y(1)-Y(3))./2.0+1.0)-(sin(Y(1)).*(4.9e+1./1.0e+1)+sin(Y(3)).*(1.47e+2./1.0e+1))./(cos(Y(1)-Y(3))./2.0+2.0)+Y(2)./(cos(Y(1)-Y(3)).*5.0+1.0e+1)+Y(4)./(cos(Y(1)-Y(3)).*5.0+2.0e+1)-((Y(1)-Y(3)).*Y(4).^2)./(cos(Y(1)-Y(3))+2.0)+(sin(Y(1)-Y(3)).*(Y(2).^2-Y(4).^2))./(cos(Y(1)-Y(3))+4.0)]

interval = [0 20];
xInit = [2 0];
yInit = [2 0];
ySol = ode45(M,interval,[xInit, yInit]);

Warning: Failure at t=2.440885e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.440892e-16) at time t.

figure
plot(ySol.x, ySol.y)
grid
legend(string(Subs), 'Location','best')

tValues = linspace(0,20,100);
yValues = deval(ySol,tValues,1);

Error using deval (line 137)
Attempting to evaluate the solution outside the interval [0.000000e+00, 2.440885e-01] where it is defined.

plot(tValues,yValues)

There are still problems, however the code essentially works.

.

[Nikhil Kapoor]

How to get values of the subs functions like y and x