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
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
James Tursa
am 9 Apr. 2021
Bearbeitet: James Tursa
am 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?
Antworten (1)
James Tursa
am 9 Apr. 2021
Bearbeitet: James Tursa
am 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 Kommentare
Siehe auch
Kategorien
Mehr zu Numerical Integration and Differential 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!