Help on solving nonlinear coupled equations with Runge-Kutta approach
20 views (last 30 days)
Show older comments
I have the following piece of code. I'm providing a "piece" here, since the full code is quite lengthy and I think the following would suffice. But I'm more than happy to share more, if needed.
x = [0, 5, 0, 0, 0, 0].';
% Simulation parameters
dt = 0.1; % time step (s)
tend = 10; % end time (s)
tspan = 0:dt:tend; % time vector
number_of_time_steps = tend/dt; % number of time steps
for t=1:number_of_time_steps
f = @(x,t) eomsys(x, data); % I have more parameters than these two.
[x, dxdt] = rungekutta4(f,x,u,dt);
end
My equations of motions are (nonlinear, coupled) quite similar to the ones I've posted in earlier posts (a quite old one here, and here), but with some additions. Basically, I have formulated the equations in following form, and the script follows. I'm sorry about the messy matrices!
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1345129/image.jpeg)
function [xdot] = eomsys(x, data)
u = x(1); v = x(2); w = x(3); p = x(4); q = x(5); r = x(6);
massmat = [70.0532949500000 0 0.715271137500000 0 2.28968497500000 0
0 93.4368680000000 0 0.445077447500000 0 -0.607912125000000
1.33294997500000 0 138.803890750000 0 0.520733722500000 0
0 0.444320895000000 0 10940.2695433200 0 0.0816402557500000
5.67678722500000 0 55153.0589243300 0 12.7759792500000 0
0 -0.608241150000000 0 0.0810982870000000 0 55162.3832866000];
smat = [0, 86.9651*r, 4.1729*w, 0.9520*r, 0, 0.9520*r;
(-2.7108*r-0*v-0*p-0*0*r), -2.3186, 0, -0.1296, 0, -11.1265;
0, 0, 0.5241, 0, -76.0312, 0;
(0*v-0*p-0*0*r), -0.3209, 0, -0.0103, 0, (1.0939e+04*q);
0, 0, 13.0970, 0, -0.0382*q, -1.0940e+04*p;
(2.7108*v - (65.4088+0*0)*v - (0.9520+0*0)*p +...
(1.1414+0*0^2)*r), -25.5156, 0, -0.6808, -4.4212e+04*p, -2.2918];
kmat = [1; 1; 1; 1; 1; 1];
% I have added ones here arbitrarily. I will have a different, and a bigger matrix here.
xdot = inv(massmat)*smat*x + inv(massmat)*kmat;
end
And my Runge-Kutta scheme,
function [xnxt, dxdt] = rungekutta4(f,x,u,dt)
xo = x;
dxdt = feval(f,xo,u);
k1 = dt*dxdt;
x = xo+0.5*k1;
k2 = dt*feval(f,x,u);
x = xo+0.5*k2;
k3 = dt*feval(f,x,u);
x = xo+k3;
k4 = dt*feval(f,x,u);
xnxt = xo + (k1+2*(k2+k3)+k4)/6;
end
As you can see, I'm learning as I go, and I feel like my approach here is wrong. I'll start from the obvious one - Can I write my nonlinear coupled ODEs like this, in matrix form? (I read here that it's not possible, but I had found a few tutorials online that did otherwise). If that is the case, this post would be redundent, and I'd appreciate a suggestion to move forward.
If I am wrong here, what am I doing wrong? If, somehow, I am right, what can I do to improve?
Thank you in advance!
11 Comments
Torsten
on 6 Apr 2023
What you solve for in the example Matlab example - Solve Equations of Motion for Baton Thrown into Air is displacement and velocity. Acceleration is what you prescribe: inv(M)*rhs .
Accepted Answer
Sam Chak
on 5 Apr 2023
Hi @Jake
I am unfamiliar with your equations of motion. But you can check and compare if you get similar plots below as solved by ode45().
% Simulation parameters
dt = 0.01; % time step (s)
tend = 10; % end time (s)
tspan = 0:dt:tend; % time vector
x0 = [0, 5, 0, 0, 0, 0];
[t, x] = ode45(@eomsys, tspan, x0);
for j = 1:6
subplot(3, 2, j)
plot(t, x(:,j)), grid on
title(sprintf('x_{%d}', j))
end
function [xdot] = eomsys(t, x)
xdot = zeros(6, 1);
u = x(1); v = x(2); w = x(3); p = x(4); q = x(5); r = x(6);
massmat = [70.0532949500000 0 0.715271137500000 0 2.28968497500000 0
0 93.4368680000000 0 0.445077447500000 0 -0.607912125000000
1.33294997500000 0 138.803890750000 0 0.520733722500000 0
0 0.444320895000000 0 10940.2695433200 0 0.0816402557500000
5.67678722500000 0 55153.0589243300 0 12.7759792500000 0
0 -0.608241150000000 0 0.0810982870000000 0 55162.3832866000];
smat = [0, 86.9651*r, 4.1729*w, 0.9520*r, 0, 0.9520*r;
(-2.7108*r-0*v-0*p-0*0*r), -2.3186, 0, -0.1296, 0, -11.1265;
0, 0, 0.5241, 0, -76.0312, 0;
(0*v-0*p-0*0*r), -0.3209, 0, -0.0103, 0, (1.0939e+04*q);
0, 0, 13.0970, 0, -0.0382*q, -1.0940e+04*p;
(2.7108*v - (65.4088+0*0)*v - (0.9520+0*0)*p +...
(1.1414+0*0^2)*r), -25.5156, 0, -0.6808, -4.4212e+04*p, -2.2918];
kmat = [1; 1; 1; 1; 1; 1];
xdot = inv(massmat)*smat*x + inv(massmat)*kmat;
end
More Answers (0)
See Also
Categories
Find more on Systems of Nonlinear Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!