Azzera filtri
Azzera filtri

Help on solving nonlinear coupled equations with Runge-Kutta approach

17 visualizzazioni (ultimi 30 giorni)
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
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 .
Jake
Jake il 7 Apr 2023
You're right, @Sam Chak :) This relates to aircraft dynamics and I'm working with a system of fully-coupeld translational and rotational motion equations.
Thank you for the response, and I understand @Torsten. Unfortunately, I think I lack the knowledge to fully comprehend all these, and I'm going to go back to learning a bit more and deriving the equations again before I give up completely.
My equations are kind of similar to the ones shown in the example (Matlab example - Solve Equations of Motion for Baton Thrown into Air), but I have a bit more going on. For example, my motion equation in x direction is in this form. The system has 6 degrees of freedom.
I have tried several approaches, but I couldn't figure a proper approach. I tried with an ode solver, but I do need intermediate values of the coefficients (such as f_1 and k_1), so I chose Runge-Kutta with a time loop. As you can guess, the kmat I have included in my above code was supposed to represent the time and velocity dependent values.
I'm sorry if I'm taking too much of time, and I can open another post if you think it's best. But I'd gladly accept any feedback at this point.

Accedi per commentare.

Risposta accettata

Sam Chak
Sam Chak il 5 Apr 2023
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)

Categorie

Scopri di più su Programming in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by