Runge-Kutta 2

1 次查看(过去 30 天)
Nagy Csaba Norbert
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 个评论
Nagy Csaba Norbert
One more question, i have a task with initialize conditions at x(0) = 1 and x'(0) = 1. Can someone also help with this?
James Tursa
James Tursa 2021-4-9
编辑:James Tursa 2021-4-9
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?

请先登录,再进行评论。

回答(1 个)

James Tursa
James Tursa 2021-4-9
编辑:James Tursa 2021-4-9
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),'.-')

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by