Solving Freedom System using Runge Kutta and Adams Bashford method

Hi I am looking to solve the above system using both Runge Kutta 4th order and Adams Bashford 4th order. The initial conditions are as follows;

x1(0)=-2

x2(0)=5

dx1/dt(0)=1

dx2/dt(0)=-3

m1=2

m2=5

k1=2

k2=4

k3=6

c1=0.06

c2=0.04

c3=0.02

Any help would be highly appreciated about where to start with the code thanks in advance.

[Sam Chak]

Hi @Lucy Martin

Good to know that you have enthusiasm about starting to write the code yourself. It will be a good learining experience.

First of all, can you provide the formulas in Runge–Kutta and Adams–Bashford 4th-order methods? I do not memorize them and it is absolutely a positive way to contribute something on your part too.

There are generally 3 parts in the numerical method:

1. the dynamical system (the governing equations derived from some fundamental laws and practical assumptions)
2. user-defined parameters (initial values, constants and customized variables that you can change)
3. the numerical solver (selection of step size, and the algorithm used to perform numerical integration of ODEs)

Item #1 is shown on the image (though you need to rearrange the ODEs for the numerical integration purposes). Item #2 is provided. For Item #3, you need to provide the step size which affects the desired accuracy of the numerical solution. The formulas for the numerical solvers must be listed out.

These are the basic preparations before writing the code.

[Lucy Martin]

function [wi, ti] = AB4 ( RHS, t0, x0, tf, N )

if size(x0,2)>1

x0=x0';

end

neqn = length ( x0 );

ti = linspace ( t0, tf, N+1 );

wi = zeros( neqn, N+1);

wi(1:neqn, 1) = x0;

h = ( tf - t0 ) / N;

oldf = zeros(neqn,N);

% Implementation of RK4

for i = 1:3

k1 = h * feval ( RHS, t0, x0 );

k2 = h * feval ( RHS, t0 + h/2, x0 + k1/2 );

k3 = h * feval ( RHS, t0 + h/2, x0 + k2/2 );

k4 = h * febal ( RHS, t0 + h, x0 + k3 );

x0 = x0 + ( k1 + 2*k2 + 2*k3 + k4 ) /6;

t0 = t0 + h;

wi(1:neqn,i+1) = x0;

end

% Implementation of Adams Bashford

for i = 4:N

oldf(1:neqn,i) = feval ( RHS, t0, x0 );

x0 = x0 + (h/24) * ( 55*oldf(:,i) - 59*oldf(:,i-1) + 37*oldf(:,i-2) - 9*oldf(:,i-3) );

t0 = t0 + h;

wi(1:neqn,i+1) = x0;

end

Thanks for the reply here is my code for the implementation of AB4 and RK4 and my desired step size is 0.25.

Answer 1

[Sam Chak]

Hi @Lucy Martin

This is the script for the mass-spring-damper system using the Runge–Kutta 4th order Solver. You can compare the result with ODE45 Solver. The script is given below as well. If you are okay, try implementing the Adams–Bashforth method yourself as a coding skill exercise. Hope this Answer is satisfactory.

function Demo_RK4

close all;

clc

global m1 m2 c1 c2 c3 k1 k2 k3

m1 = 2;

m2 = 5;

c1 = 0.06;

c2 = 0.04;

c3 = 0.02;

k1 = 2;

k2 = 4;

k3 = 6;

tStart = 0;

step = 0.25; % if step is smaller, then numerical solution is more accurate

tEnd = 20;

t = [tStart:step:tEnd]; % simulation time span

y0 = [-2 5 1 -3]; % initial values

f = @(t,x) [x(3); % system of 1st-order ODEs

x(4);

- ( (k1+k2)*x(1) - k2*x(2) + (c1+c2)*x(3) - c2*x(4))/m1;

- (- k2*x(1) + (k2+k3)*x(2) - c2*x(3) + (c2+c3)*x(4))/m2];

yRK4 = RK4Solver(f, t, y0); % calling Runge-Kutta 4th-order Solver

% plotting the numerical solution

plot(t, yRK4, 'linewidth', 1.5)

grid on

xlabel('Time, t [sec]')

ylabel({'$x_{1},\; x_{2},\; \dot{x}_{1},\; \dot{x}_{2}$'}, 'Interpreter', 'latex')

title('Time responses of the system states (with RK4 Solver)')

legend({'$x_{1}$', '$x_{2}$', '$\dot{x}_{1}$', '$\dot{x}_{2}$'}, 'Interpreter', 'latex', 'location', 'best')

end

function y = RK4Solver(f, x, y0)

y(:, 1) = y0; % initial condition

h = x(2) - x(1); % step size

n = length(x); % number of steps

for j = 1 : n-1

k1 = f(x(j), y(:, j)) ;

k2 = f(x(j) + h/2, y(:, j) + h/2*k1) ;

k3 = f(x(j) + h/2, y(:, j) + h/2*k2) ;

k4 = f(x(j) + h, y(:, j) + h*k3) ;

y(:, j+1) = y(:, j) + h/6.0*(k1 + 2*k2 + 2*k3 + k4) ;

end

end

Result:

Here is the script for the mass-spring-damper system using the ODE45 Solver.

function Demo_ode45

close all

clc

tspan = [0 20]; % time span of simulation

y0 = [-2 5 1 -3]; % initial values

[t, y] = ode45(@(t, y) odefcn(t, y), tspan, y0);

plot(t, y, 'linewidth', 1.5)

grid on

xlabel('Time, t [sec]')

ylabel({'$x_{1},\; x_{2},\; \dot{x}_{1},\; \dot{x}_{2}$'}, 'Interpreter', 'latex')

title('Time responses of the system states (with ODE45 Solver)')

legend({'$x_{1}$', '$x_{2}$', '$\dot{x}_{1}$', '$\dot{x}_{2}$'}, 'Interpreter', 'latex', 'location', 'best')

end

function dydt = odefcn(t, y) % Ordinary Differenctial Equations

dydt = zeros(4,1);

m1 = 2;

m2 = 5;

c1 = 0.06;

c2 = 0.04;

c3 = 0.02;

k1 = 2;

k2 = 4;

k3 = 6;

M = diag([m1 m2]);

C = [c1+c2 -c2; -c2 c2+c3];

K = [k1+k2 -k2; -k2 k2+k3];

dydt(1) = y(3);

dydt(2) = y(4);

dydt(3:4) = M\(- K*[y(1); y(2)] - C*[y(3); y(4)]);

end

Result:

[Lucy Martin]

Thanks a lot this code worked perfectly have now implemented the AB4 method and it works great too.

[Sam Chak]

Hi @Lucy Martin

I'm glad that your AB4 works great! If you like my Answer, please Vote and Accept. 😄