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

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

[Torsten]

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]

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)

Answer 1

[Nipun]

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