Help on solving nonlinear coupled equations with Runge-Kutta approach
20 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Torsten
am 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 .
Akzeptierte Antwort
Sam Chak
am 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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Systems of Nonlinear Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!