Need assistance with incorrect output in plot.

1 view (last 30 days)
Hi all. I am working on a class project that takes a set of four first order equations and integrates them into a plot. The plot should appear as this below but my graph is totally different and I am not sure what is causing the issue. Can anyone point out where I am going wrong? I am including my full code below. Thanks in advance
%% Creating custom R-K Integrator
clear
clc
dt=0.0001 %0.01 sec was used initially,
x_vect=[0,0,0,0]; %Initial Conditions
t_span=[0:dt:2];
npoints=length(t_span);
x_dot = @deriv
t0=0;
t(1)=t0;
x_answer(1)=x_vect(1,4);% Initialize answer vector.
for i=1:npoints
t(i+1) = t(i)+dt; %update t, this value will not be used until next iteration.
x_vect = integ(x_vect,x_dot,t(i),dt);
x_answer(i+1)=x_vect(1,4);
t0=t(i);
end
plot(t,x_answer)
function out=integ(x_vect,x_dot,t,dt) %fourth order runge kutta
%out=x_vect+x_dot*dt;
k1=x_dot(t,x_vect)*dt;
k2=x_dot(t+dt/2,x_vect+k1/2)*dt;
k3=x_dot(t+dt/2,x_vect+k2/2)*dt;
k4=x_dot(t+dt,x_vect+k3)*dt;
out=x_vect+1/6*(k1+2*k2+2*k3+k4);
end
function out=deriv(t,y_init)
out(1) = y_init(2);
out(2) = y_init(3);
out(3) = y_init(4);
out(4) = -206*y_init(4)-1341.4*y_init(3)-848.4*y_init(2)+60000;
end
  5 Comments
darova
darova on 29 Mar 2021
x_dot according to my reacher is supposed to represent the derivitive.

Sign in to comment.

Answers (1)

darova
darova on 29 Mar 2021
Try to change this part
x_answer = zeros(npoints,4); % preallocate answer vector
x_answer(1,:) = x_vect(1,4); % Initialize answer vector.
And this part inside for loop
x_answer(i+1,:) = x_vect
  1 Comment
Sean Sweeney
Sean Sweeney on 29 Mar 2021
x_dot = @deriv
t0=0;
t(1)=t0;
x_answer = zeros(npoints,4);
x_answer(1,:) = x_vect(1,4);
for i=1:npoints
t(i+1) = t(i)+dt; %update t, this value will not be used until next iteration.
x_vect = integ(x_vect,x_dot,t(i),dt);
x_answer(i+1,:) = x_vect
t0=t(i);
end
The bold is what I updated from the original code, but it just ran forever until I manually ended the script. Was this the change you were suggesting?

Sign in to comment.

Categories

Find more on Programming 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!