Solving a System of 2nd Order Nonlinear ODEs

Hi,

I am having some issues plotting this two equations. For some reason MATLAB is taking forever to solve. Is my syntax incorrect? See attached equation

function dydt = odefcn(t, y)    
% y1 = x
% y2 = theta
% y3 = xdot
% y4 = thetadot
    dydt = zeros(4,1);
    dydt(1) = y(3); % linear speed
    dydt(2) = y(4); % angular speed
    dydt(3) = ((6*sin(y(2)))/((-3*(sin(y(2)))^2)+8))*((-1.5/2)*sin(y(2)) - (4/(3*sin(y(2))))*(-(cos(y(2))/4)*(y(4))^2 + 0.1*y(3) + y(1))); % linear acceleration
    dydt(4) = ((((-6*sin(y(2)))/(-3*(sin(y(2)))^2))*((-1.5/2)*sin(y(2)) - (4/(3*sin(y(2))))*(-(cos(y(2))/4)*(y(4))^2 + 0.1*y(3) + y(1)))) - ((cos(y(2)))/4)*(y(4))^2 + 0.1*y(3) + y(1))*(4/sin(y(2))); % angular acceleration
end

function mpz
close all
clc  
tspan = [0 2];                 % time span of simulation
y0 = [0.5 pi/3 0 0];            % initial values
[t, y] = ode45(@(t, y) odefcn(t, y), tspan, y0);
plot(t, y(1), 'linewidth', 1.5)
grid on
xlabel('Time, t [sec]')
ylabel({'$x,\; \theta,\; \dot{x},\; \dot{\theta}$'}, 'Interpreter', 'latex')
title('Time responses of the system states')
legend({'x', '$\theta$', '$\dot{x}$', '$\dot{\theta}$'}, 'Interpreter', 'latex', 'location', 'best')
end

Answer 1

[Sam Chak]

Hi @mpz

A few days ago, you have learned how to uncouple the ODEs by using the substitution method in this link:

https://www.mathworks.com/matlabcentral/answers/1687154-solving-a-system-of-2nd-order-nonlinear-odes

This time, you will learn the matrix method that is conceptually similar to the inverse method described by @William Rose, but the matrix size is 75% smaller. Computing the determinant for a matrix is also much simpler than a matrix.

Let's begin! Equations 7 and 8 can be rewritten in matrix form

.

Solving for and

.

Now you can write the ODEs separately as:

function mpz

close all

clc

tspan = [0 40]; % time span of simulation

y0 = [0.5 pi/3 0 0]; % initial values

[t, y] = ode45(@(t, y) odefcn(t, y), tspan, y0);

plot(t, y(:,1), t, y(:,2), 'linewidth', 1.5)

hold on

plot(t, y(:,3), t, y(:,4))

hold off

grid on

xlabel('Time, t [sec]')

ylabel({'$x,\; \theta,\; \dot{x},\; \dot{\theta}$'}, 'Interpreter', 'latex')

title('Time responses of the system states')

legend({'x', '$\theta$', '$\dot{x}$', '$\dot{\theta}$'}, 'Interpreter', 'latex', 'location', 'best')

end

function dydt = odefcn(t, y) % Ordinary Differential Equations

dydt = zeros(4,1);

dydt(1) = y(3);

dydt(2) = y(4);

dydt(3) = (1/((1/8)*(sin(y(2)))^2 - 1/3))*((3/16)*(sin(y(2)))^2 + (1/3)*(y(1) + (1/10)*y(3) - (1/4)*cos(y(2))*(y(4))^2));

dydt(4) = (1/((1/8)*(sin(y(2)))^2 - 1/3))*((3/4)*sin(y(2)) + (1/2)*sin(y(2))*(y(1) + (1/10)*y(3) - (1/4)*cos(y(2))*(y(4))^2));

end

Results:

From the plot, we can see that x is converging faster than θ.

[William Rose]

@Sam Chak, that was very nice!

[Sam Chak]

Thank you, @William Rose.

I actually learned something new about the ode45 capability from the link you posted:

https://www.mathworks.com/help/matlab/math/solve-equations-of-motion-for-baton.html

It is very interesting because it saves time for some people (doing the matrix inverse is actually quite tedious for coupled systems, especially in some motion systems with 6 and higher-degree-of-freedom).

Also, I wonder if it works for the Descriptor Systems like:

with can possibly become singular, as well as with a regular matrix pencil.

[mpz]

Thank you @William Rose and @Sam Chak

Answer 2

[William Rose]

@mpz,

Please show how you converted the two differential equations, labeleled "7" and "8" in your figure, to a set of first order O.D.E.s. I think that this set of ODEs is implicit, i.e. you will have to use a mass matrix in the solution. You can still use ode45() with an implicit set of equations. Specify the mass matrix using the Mass option of odeset. See the help for details.

[mpz]

@William Rose

see attached on how I got it. I used non-dimensionalization to convert. I have looked at it multiple times but can't see what is incorrect. I haven't used the odeset before.

[mpz]

@William Rose

this is the set

[William Rose]

[edit: fix error in eq. 10, eq.12, M]

Define:

Then equation 7 becomes

(Eq.9)

Equation 8 becomes

(Eq.10)

Rewrite 9:

(Eq.11)

Rewrite 10:

(Eq.12)

Now we have

where and f(y) are column vectors, and

,

See https://www.mathworks.com/help/matlab/math/solve-equations-of-motion-for-baton.html for how to code this. Figure out the equations yourself and see if you agree with my results above.

[William Rose]

@mpz,

Maybe you have already inverted the mass matrix and multiplied . If so, your way should work, and you don't have to use the mass matix option.

I see I made a mistake in equation 10 above. That mistake propagates into equation 12 and the equation for matrix M.

[mpz]

@William Rose I was coded the equations but getting a failure.

function M = mass(t,y)
% Mass matrix function
M = zeros(4,4);
M(1,1) = 1;
M(2,2) = 1;
M(3,3) = -0.5*sin(y(2));
M(3,4) = 1/3;
M(4,3) = 1;
M(4,4) = -sin(y(2));
end

function dydt = f(t,y)
% Equation to solve
dydt = [y(3)
        y(4)
        -0.75*sin(y(2))
        0.25*cos(y(2))*(y(4))^2 - 0.1*y(3) - y(1)];
end

tspan = linspace(0,20,1000);
y0 = [0.5 pi/3 0 0];
opts = odeset('Mass',@(t,y) mass(t,y));
[t,y] = ode45(@(t,y) f(t,y),tspan,y0,opts);

Here is the stop warning text

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

[William Rose]

My mass matrix equation had an error in M(4,4), which I have now fixed in my comment above. I omitted a factor of 0.25 from M(4,4). If you update your code to include the corrected M(4,4), does it run better?

Try setting the Mass state dependence to 'strong' - see help for odeset. Example:

opts = odeset('Mass',@M,'MStateDependence','strong');

If you still get errors or have problems, you could try a different ODE solver. Review the recommendations here. Maybe ode15s or ode15i.