One more question, i have a task with initialize conditions at x(0) = 1 and x'(0) = 1. Can someone also help with this?
Runge-Kutta 2
2 views (last 30 days)
Show older comments
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
on 9 Apr 2021
Edited: 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?
Answers (1)
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),'.-')
0 Comments
See Also
Categories
Find more on Numerical Integration and Differential Equations 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!