Help on solving nonlinear coupled equations with Runge-Kutta approach
17 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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!
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 Commenti
Torsten
il 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 .
Risposta accettata
Sam Chak
il 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
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!