Solving a System of 2nd Order Nonlinear ODEs
조회 수: 25 (최근 30일)
이전 댓글 표시
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
댓글 수: 0
채택된 답변
Sam Chak
2022년 4월 8일
편집: Sam Chak
2022년 4월 8일
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 θ.
댓글 수: 3
Sam Chak
2022년 4월 8일
I actually learned something new about the ode45 capability from the link you posted:
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.
추가 답변 (1개)
William Rose
2022년 4월 8일
@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.
댓글 수: 6
William Rose
2022년 4월 8일
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.
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!