MATLAB Answers

Runge-Kutta 2

4 views (last 30 days)
Nagy Csaba Norbert
Nagy Csaba Norbert on 7 Apr 2021
Edited: James Tursa on 9 Apr 2021
So I'm doing a project for school, where i have to solve the x"+5x'+4x = 3 - 2t - t^2 equation with Runge-Kutta 2.
intervalmin = 0;
intervalmax = 1;
h = 0.1;
g = 0.02;
numnodes = ceil(((intervalmax - intervalmin) / h) + 1);
numnodes2 = ceil(((intervalmax - intervalmin) / g) + 1);
inival = 0;
t = zeros(1, numnodes);
x = zeros(2, numnodes);
t(1) = intervalmin;
x(:,1) = inival;
f = @(t,x) [x(2); -5 * x(2) - 4 * x(1) + 3 - 2 * (t) - t.^2];
for i = 2:numnodes
t(i) = t(i - 1) + h;
k1 = f(t(i - 1), x(:, i - 1));
k2 = f(t(i - 1) + h / 2, x(:, i-1) + (h / 2) *k1);
k3 = f(t(i - 1) + h, x(:, i - 1) - h * k1 + 2 * h * k2);
x(:, i) = x(:, i - 1) + (h / 6) * (k1 + 4 * k2 + k3);
end
o = zeros(1, numnodes2);
z = zeros(2, numnodes2);
o(1) = intervalmin;
z(:, 1) = inival;
m = @(o, z) [z(2); -5 * z(2) - 4 * z(1) + 3 - 2 * (o) - o.^2];
for i = 2:numnodes2
o(i) = o(i - 1) + g;
k1 = m(o(i - 1), z(:, i - 1));
k2 = m(o(i - 1) + h / 2, z(:, i - 1) + (g / 2) * k1);
k3 = m(o(i - 1) + h, z(:, i - 1) - g * k1 + 2 * g * k2);
z(:, i)=z(:, i - 1) + (g / 6) * (k1 + 4 * k2 + k3);
end
[TT, XX] = ode45(f,[intervalmin intervalmax],x(:, 1));
figure
plot(t,x(2,:),'.-',o,z(2,:),'.-',TT,XX(:,2),'.-')
hold on
syms x(t)
dz = diff(x);
ode = diff(x,t,2) == (-5 * x(2) - 4 * x(1) + 3 - 2 * (t) - t.^2);
cond1 = x(0) == 1;
cond2 = dz(0) == 1;
dsolve(ode, cond1, cond2)
legend('h=0.1','h=0.02','ode45')
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
grid on;
I'm pretty new to matlab, and just can't seem to find a way to make it work, so if somebody could help me solve it, i would be grateful
  2 Comments
James Tursa
James Tursa on 9 Apr 2021
This looks like a 3rd order method your are implementing. Is that what you are supposed to be doing? Or are you supposed to be implementing a 2nd order method as your title suggests?

Sign in to comment.

Answers (1)

James Tursa
James Tursa on 9 Apr 2021
Edited: James Tursa on 9 Apr 2021
Normally one would plot the position, not the velocity. So I would have expected you to plot the "1" index of your solutions. E.g.,
plot(t,x(1,:),'.-',o,z(1,:),'.-',TT,XX(:,1),'.-')

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by