How to implement Forward Euler Method on a system of ODEs

Hi,

I was wondering if it is possible to solve a function using the forward Euler method in the form:

function [dydt] = ODEexample(t,y, AdditionalParameters)

dydt(1:2) = y(3:4);
dydt(3:4) = y(1:2)
dydt = dydt';
end

Using a Euler method such as https://nl.mathworks.com/matlabcentral/answers/366717-implementing-forward-euler-method, .

Or should I adjust the function?

[darova]

Yes, it's possible. Do you have any attempts?

[Cynthia Struijk]

Yes I do. I created a function which is kinda analogous to the link I mentioned, and it already get an error for the fourth line for reshaping yold, the inital condition:

function [tnew, ynew] = FirstOrderEulerSystem(fun, tspan, yold, par ,n)
% Compute the Ordinary differential Equation solution vectors and
% corresponding time
% Inputs: (function, time span, initial values for the ODE, additional
% parameters, time steps)
% Outputs: time and solution ODE (as time steps x 2 matrix)
% Example: 
%[tEuler, yEuler] = FirstOrderEulerSystem(@ODEsystem, [0 2000], [12 344], par, 100);
  %The number of equations can be determined by using the length command
  %over the initial conditions (at least one initial condition is given for each
  %equation)
  NumberOfEquations = length(yold);
  
  %Determine the step size h, by dividing the difference between the final
  %and starting time by the amount of the time steps n
  h = (tspan(2) - tspan(1))/n;

  %First values for t and y 
  told = tspan(1);  %initial time
  yold = reshape(yold, 1, NumberOfEquations); %initial condition, reshaped in 1xNumberofEquation matrix
  

When I tried to implement a function like ODEexample in this function like:

%[tEuler, yEuler] = FirstOrderEulerSystem(@ODEexample, [0 2000], [ [1;0] [0;-4.5] ], par, 100);

I get the error:

Error using reshape
To RESHAPE the number of elements must not change.

Error in FirstOrderEulerSystem (line 21)
  yold = reshape(yold, 1, NumberOfEquations); %initial condition, reshaped in 1xNumberofEquation matrix
  

I guess it has something to do with the fact that the function ODEexample and the initial conditions are in a form of a matrix, but I am not sure if this is true?

[darova]

Where are the original equations?

[Cynthia Struijk]

 function [dxdt] = ODEparticles(t,x, par_c)
dxdt(1:2) = x(3:4);
dxdt(3:4) = (par_c.k*par_c.qi*par_c.qj*(x(1:2)-par_c.xj))/(par_c.m*(norm(par_c.xj-x(1:2))^3));
dxdt = dxdt';
 end

with

 par_c.k = 14.39964548; % electronvolt per Angström
par_c.xj = [0;0]; % position of particle j
par_c.m = 1; % mass of particle in unified atomic mass units  
par_c.qi = 1; % charge of electron i in Coulomb
par_c.qj = -1; % charge of electron j in Coulomb

tspan_particle = [0 20];
init_particle = [ [1;0] [0;-4.5] ];

Answer 1

[darova]

I've tried and reached a success

par_c.k = 14.39964548; % electronvolt per Angström
par_c.xj = [0 0]; % position of particle j
par_c.m = 1; % mass of particle in unified atomic mass units  
par_c.qi = 1; % charge of electron i in Coulomb
par_c.qj = -1; % charge of electron j in Coulomb

T = 2;                  % simulation time
n = 100;                % number of steps
X = zeros(n,4);         % preallocation
dt = T/n;               % time step
X(1,:) = [ 1 0 0 -4.5 ];% initial conditions?
 
f = @(x) [x(3:4) (par_c.k*par_c.qi*par_c.qj*(x(1:2)-par_c.xj))/(par_c.m*(norm(par_c.xj-x(1:2))^3))];
for i = 1:n-1
    X(i+1,:) = X(i,:) + dt*f(X(i,:));
end
plot(X)

[Cynthia Struijk]

The theoretical solution and the ode45 solution predict a circular motion. Using this code below:

tspan_particle = [0 20];
init_particle = [ 1 0 0 -4.5 ];
options_c = odeset('RelTol', 1e-4, 'AbsTol', [1e-4 1e-4 1e-4 1e-4]);
[t_particle,x_particle] = ode45(@ODEparticles ,tspan_particle ,init_particle ,options_c, par_c);

The Forward Euler method has to give the same result right?

[darova]

Try to increase number of steps

n = 1000;                % number of steps

[Cynthia Struijk]

Unfortunalely, even a million steps would not work :( I really appreciate your help though!

[darova]

Impossible!

function main

par_c.k = 14.39964548; % electronvolt per Angstr?m
par_c.xj = [0 0]; % position of particle j
par_c.m = 1; % mass of particle in unified atomic mass units
par_c.qi = 1; % charge of electron i in Coulomb
par_c.qj = -1; % charge of electron j in Coulomb

T = 5;                  % simulation time
n = 1000;                % number of steps
X = zeros(n,4);         % preallocation
dt = T/n;               % time step
X(1,:) = [ 1 0 0 -4.5 ];% initial conditions?

    function dy = f(~,x)
        x = x(:)';
        dy = [x(3:4) (par_c.k*par_c.qi*par_c.qj*(x(1:2)-par_c.xj))/(par_c.m*(norm(par_c.xj-x(1:2))^3))];
        dy = dy';
    end
for i = 1:n-1
    X(i+1,:) = X(i,:) + dt*f(1,X(i,:))';
end
plot(X(:,1),X(:,2))

[t,X] = ode45(@f,[0 T],[1 0 0 -4.5]);
hold on
plot(X(:,1),X(:,2),'r')
hold off
legend('euler method','ode45')
end

[Cynthia Struijk]

Oh probably I did something wrong then. This seems to work, thank you so much for all the help, hope you'll have a wonderful day :)