How can I fix my Runge Kutta (4th order) method to solve a 2nd order ODE?

5 次查看（过去 30 天）
Mert
Mert 2024-5-20
回答： Nipun 2024-6-3
Ran in:
solve x''(t) +δx'(t) + αx(t) + βx(t)^3 = γcos(ωt), x(0)=0 , x'(0)=0
clc
clear;
h=1;
delta = 0.1;
alpha = -1;
beta = 0.25;
gamma = 2.5;
omega = 2;
t(1)=0;
z(1)=0;
x(1)=0;
d=(gamma*(cos(omega*t)))-(beta*(x^3))-(alpha*x)-(delta*z)
d = 2.5000
g=@(t,x,z) z;
f=@(t,x,z) (gamma*(cos(omega*t)))-(beta*(x^3))-(alpha*x)-(delta*z);
for i=1:50
L1 = h*g(t, x, z);
k1 = h*f(t, x, z);
L2 = h*g(t+h/2, x+k1/2, z+L1/2);
k2 = h*f(t+h/2, x+k1/2, z+L1/2);
L3 = h*g(t+h/2, x+k2/2, z+L2/2);
k3 = h*f(t+h/2, x+k2/2, z+L2/2);
L4 = h*g(t+h, x+k3, z+L3);
k4 = h*f(t+h, x+k3, z+L3);
z = z + (L1+2*L2+2*L3+L4)/6;
x = x + (L1+2*L2+2*L3+L4)/6;
t= t+h;
fprintf('i=%8.0f, t=%8.2f, x=%8.6f, z=%8.6f\n',i,t,x,z)
plot(t,z,'-*r')
hold on
plot(t,x,'-ob')
end
i= 1, t= 1.00, x=0.000000, z=0.000000 i= 2, t= 2.00, x=0.000000, z=0.000000 i= 3, t= 3.00, x=0.000000, z=0.000000 i= 4, t= 4.00, x=0.000000, z=0.000000 i= 5, t= 5.00, x=0.000000, z=0.000000 i= 6, t= 6.00, x=0.000000, z=0.000000 i= 7, t= 7.00, x=0.000000, z=0.000000 i= 8, t= 8.00, x=0.000000, z=0.000000 i= 9, t= 9.00, x=0.000000, z=0.000000 i= 10, t= 10.00, x=0.000000, z=0.000000 i= 11, t= 11.00, x=0.000000, z=0.000000 i= 12, t= 12.00, x=0.000000, z=0.000000 i= 13, t= 13.00, x=0.000000, z=0.000000 i= 14, t= 14.00, x=0.000000, z=0.000000 i= 15, t= 15.00, x=0.000000, z=0.000000 i= 16, t= 16.00, x=0.000000, z=0.000000 i= 17, t= 17.00, x=0.000000, z=0.000000 i= 18, t= 18.00, x=0.000000, z=0.000000 i= 19, t= 19.00, x=0.000000, z=0.000000 i= 20, t= 20.00, x=0.000000, z=0.000000 i= 21, t= 21.00, x=0.000000, z=0.000000 i= 22, t= 22.00, x=0.000000, z=0.000000 i= 23, t= 23.00, x=0.000000, z=0.000000 i= 24, t= 24.00, x=0.000000, z=0.000000 i= 25, t= 25.00, x=0.000000, z=0.000000 i= 26, t= 26.00, x=0.000000, z=0.000000 i= 27, t= 27.00, x=0.000000, z=0.000000 i= 28, t= 28.00, x=0.000000, z=0.000000 i= 29, t= 29.00, x=0.000000, z=0.000000 i= 30, t= 30.00, x=0.000000, z=0.000000 i= 31, t= 31.00, x=0.000000, z=0.000000 i= 32, t= 32.00, x=0.000000, z=0.000000 i= 33, t= 33.00, x=0.000000, z=0.000000 i= 34, t= 34.00, x=0.000000, z=0.000000 i= 35, t= 35.00, x=0.000000, z=0.000000 i= 36, t= 36.00, x=0.000000, z=0.000000 i= 37, t= 37.00, x=0.000000, z=0.000000 i= 38, t= 38.00, x=0.000000, z=0.000000 i= 39, t= 39.00, x=0.000000, z=0.000000 i= 40, t= 40.00, x=0.000000, z=0.000000 i= 41, t= 41.00, x=0.000000, z=0.000000 i= 42, t= 42.00, x=0.000000, z=0.000000 i= 43, t= 43.00, x=0.000000, z=0.000000 i= 44, t= 44.00, x=0.000000, z=0.000000 i= 45, t= 45.00, x=0.000000, z=0.000000 i= 46, t= 46.00, x=0.000000, z=0.000000 i= 47, t= 47.00, x=0.000000, z=0.000000 i= 48, t= 48.00, x=0.000000, z=0.000000 i= 49, t= 49.00, x=0.000000, z=0.000000 i= 50, t= 50.00, x=0.000000, z=0.000000
  2 个评论
Torsten
Torsten 2024-5-20
Ran in:
For comparison:
h=0.01;
delta = 0.1;
alpha = -1;
beta = 0.25;
gamma = 2.5;
omega = 2;
fun = @(t,y)[y(2);-delta*y(2)-alpha*y(1)-beta*y(1)^3+gamma*cos(omega*t)];
tspan = 0:h:5000*h;
y0 = [0;0];
[T,Y] = ode45(fun,tspan,y0);
plot(T,Y)
Sam Chak
Sam Chak 2024-5-20
Hi @Mert
Could you please clarify how you would like the for-loop to function? Specifically, where would you like the iteration variable 'i' to be inserted?
for i=1:50
L1 = h*g(t, x, z);
k1 = h*f(t, x, z);
L2 = h*g(t+h/2, x+k1/2, z+L1/2);
k2 = h*f(t+h/2, x+k1/2, z+L1/2);
L3 = h*g(t+h/2, x+k2/2, z+L2/2);
k3 = h*f(t+h/2, x+k2/2, z+L2/2);
L4 = h*g(t+h, x+k3, z+L3);
k4 = h*f(t+h, x+k3, z+L3);
z = z + (L1+2*L2+2*L3+L4)/6;
x = x + (L1+2*L2+2*L3+L4)/6;
t= t+h;
fprintf('i=%8.0f, t=%8.2f, x=%8.6f, z=%8.6f\n',i,t,x,z)

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Nipun
Nipun 2024-6-3
Hi Mert,
I understand that you want to solve the differential equation "x''(t) + δx'(t) + αx(t) + βx(t)^3 = γcos(ωt)" with initial conditions "x(0) = 0" and "x'(0) = 0".
Here are some assumptions and corrections for your MATLAB code:
  1. I assume that you are using the Runge-Kutta method for solving the differential equation.
  2. The differential equations are coupled: "x'(t) = z(t)" and "z'(t) = γcos(ωt) - βx(t)^3 - αx(t) - δz(t)".
  3. The code needs to correctly update the time "t", the state variables "x" and "z", and use arrays to store the results for plotting.
Here's a corrected and complete version of your MATLAB code:
clc;
clear;
% Parameters
h = 1;
delta = 0.1;
alpha = -1;
beta = 0.25;
gamma = 2.5;
omega = 2;
% Initial conditions
t(1) = 0;
z(1) = 0;
x(1) = 0;
% Define the functions for the differential equations
g = @(t, x, z) z;
f = @(t, x, z) (gamma * cos(omega * t)) - (beta * x^3) - (alpha * x) - (delta * z);
% Time-stepping using the 4th-order Runge-Kutta method
for i = 1:50
L1 = h * g(t(i), x(i), z(i));
k1 = h * f(t(i), x(i), z(i));
L2 = h * g(t(i) + h/2, x(i) + L1/2, z(i) + k1/2);
k2 = h * f(t(i) + h/2, x(i) + L1/2, z(i) + k1/2);
L3 = h * g(t(i) + h/2, x(i) + L2/2, z(i) + k2/2);
k3 = h * f(t(i) + h/2, x(i) + L2/2, z(i) + k2/2);
L4 = h * g(t(i) + h, x(i) + L3, z(i) + k3);
k4 = h * f(t(i) + h, x(i) + L3, z(i) + k3);
x(i+1) = x(i) + (L1 + 2*L2 + 2*L3 + L4) / 6;
z(i+1) = z(i) + (k1 + 2*k2 + 2*k3 + k4) / 6;
t(i+1) = t(i) + h;
end
% Plotting the results
figure;
plot(t, x, '-ob', 'DisplayName', 'x(t)');
hold on;
plot(t, z, '-*r', 'DisplayName', 'z(t)');
xlabel('Time t');
ylabel('State variables');
legend;
title('Solution of the Differential Equation');
grid on;
This code correctly implements the Runge-Kutta method for the given differential equations and plots the results. The "x" and "z" arrays store the state variables, and the "t" array stores the time steps.
Hope this helps.
Regards,
Nipun

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by