Solving a System of 2nd Order Nonlinear ODEs

5 views (last 30 days)
mpz
mpz on 8 Apr 2022
Commented: mpz on 10 Apr 2022
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

Accepted Answer

Sam Chak
Sam Chak on 8 Apr 2022
Edited: Sam Chak on 8 Apr 2022
Hi @mpz
A few days ago, you have learned how to uncouple the ODEs by using the substitution method in this link:
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 θ.
  3 Comments

Sign in to comment.

More Answers (1)

William Rose
William Rose on 8 Apr 2022
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.
  6 Comments
William Rose
William Rose on 8 Apr 2022
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by