Using Matlab Runge Kutta Routine
39 views (last 30 days)
Show older comments
I'm having trouble with the following code. I'm trying to test it to get it to produce the final value as it proceeds through the count and increments the count, but it isn't working. It won't go past the first increment. I'm new to Matlab (once again). Can anyone assist me?
handle=[];
%Initial Conditions
pn0 = 0.0; %Initial north position
pe0 = 0.0; %initial east position
pd0 = -100.0; %initial east position
u0 =25.0; %initial velocity along body x-axis
v0 = 0.0; %initial velocity along body y-axis
w0 = 0.0; %initial velocity along body z-axis
phi0 = 0.0; %initial roll angle
theta0 = 0.0; %initial pitch angle
psi0 = 0.0; %initial yaw rate
l0 = 0; %initial moment about ib axis
m0 = 0; %initial moment about jb axis
n0 = 0; %initial moment about kb axis
p0 = 0; %initial roll rate along ib in Fb
q0 = 0; %initial pitch rate along jb in Fb
r0 = 0; %initial yaw rate along kb in Fb
Va0 = (u0^2 + v0^2 + w0^2)^0.5; %initial velocity
%Initial forces producing acceleration - fx, fy, fz
fx = 0; %Force in the ib direction.
fy = 0; %Force in the jb direciton.
fz = 0; %Force in the kb direction.
%Plane is initially aligned so that ib,jb, and kb
%correspond to x, y, z.
%Physical Parameters
mass = 11; % kg
Jx = 0.8244; %kg m^2
Jy = 1.135;
Jz = 1.759;
Jxz = 0.1204;
S_wing = 0.55;
b = 2.8956;
c = 0.18994;
Sprop = 0.2027;
rho = 1.2682;
e = 0.9;
AR = (b^2)/S_wing;
gravity = 9.81;
%Gamma Values
G = Jx * Jz - (Jxz^2);
G1 = (Jxz * (Jx - Jy + Jz))/ G ;
G2 = (Jz * (Jz - Jy) + (Jxz^2))/ G;
G3 = Jz / G;
G4 = Jxz / G;
G5 = (Jz - Jx) / Jy;
G6 = Jxz / Jy;
G7 = ((Jx - Jy) * Jx + (Jxz ^ 2)) / G;
G8 = Jx / G;
%Establish the parameters to pass to the subroutines
p = p0;
q = q0;
r = r0;
l = l0;
m = m0;
n = n0;
u = u0;
v = v0;
w = w0;
phi = phi0;
theta = theta0;
psi = psi0;
dt = 0.1;
fxone = fx;
fyone = fy;
fzone = fz;
fxtwo = 0;
fytwo = 0;
fztwo = 0;
n=0;
numsteps=100;
total_time = 100; % total simulation time in seconds
time_step = 0.1; % time step in seconds
num_steps = total_time / time_step; % number of steps
for count = 1:numsteps
clc;
clear all;
%Runge Kutta 4 attempt;
%fx = 5*t/mass
r=0
v=0;
q=0;
w=0;
mass = 10;
tstart = (count-1)*time_step;
tstop = tstart + 0.01;
tspan = [tstart,0.01,tstop];
y0=0.0;
[t,y]=ode45(@(t,y) (r*v-q*w)+5*t/mass, tspan, y0);
plot(t,y,'-o')
formatSpec = 'The next element of the Y array is equal to %6.4f at %6.4.';
fprintf (formatSpec,y(end,:),t(end,:));
pause;
end
1 Comment
James Tursa
on 11 Aug 2024 at 0:03
You really intend to loop through 100 plots, pausing between each one?
Answers (1)
John D'Errico
on 9 Aug 2024 at 20:32
Edited: John D'Errico
on 9 Aug 2024 at 20:35
"IT WON'T GO PAST THE FIRST INCREMENT".
Lol. Exactly what do you think this does in the very beginning of the loop?
for count = 1:numsteps
clc;
clear all;
...
You spent all that time setting up all of those constants and parameters. Then the very first thing you do is delete them all.
Any guesses why it won't go any further? Whether your code is correct besides that, I have not looked. But deleting everything you did is unlikely to be a good idea. It could have been worse I suppose. You could have put an exit command in your loop. ;-)
3 Comments
John D'Errico
on 9 Aug 2024 at 21:30
Edited: John D'Errico
on 9 Aug 2024 at 21:31
In case it is? What do you think "clear all" does?
CLEAR ALL removes all variables, globals, functions and MEX links.
How can your code proceed past that point?
help clear
Walter Roberson
on 9 Aug 2024 at 21:38
Oddly enough, if you are within a for loop, "clear all" erases the loop control variable for the rest of the duration of that iteration, but the for loop recreates the variable the next iteration.
for k = 1 : 3; if k == 1; clear all; else; disp(k); end; end
for k = 1 : 3; if k == 1; clear all; end; disp(k); end
so the variable was recreated for iterations 2 and 3.
See Also
Categories
Find more on Introduction to Installation and Licensing 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!