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

1 Ansicht (letzte 30 Tage)
Mert
Mert am 20 Mai 2024
Beantwortet: Nipun am 3 Jun. 2024
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 Kommentare
Torsten
Torsten am 20 Mai 2024
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 am 20 Mai 2024
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)

Melden Sie sich an, um zu kommentieren.

Antworten (1)

Nipun
Nipun am 3 Jun. 2024
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

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Produkte


Version

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by