MATLAB Answers

How to implement Forward Euler Method on a system of ODEs

104 views (last 30 days)
Cynthia Struijk
Cynthia Struijk on 13 Dec 2019
Edited: Cynthia Struijk on 14 Dec 2019
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
Or should I adjust the function?

  4 Comments

Show 1 older comment
Cynthia Struijk
Cynthia Struijk on 14 Dec 2019
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?
Cynthia Struijk
Cynthia Struijk on 14 Dec 2019
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] ];

Sign in to comment.

Accepted Answer

darova
darova on 14 Dec 2019
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)

  5 Comments

Show 2 older comments
Cynthia Struijk
Cynthia Struijk on 14 Dec 2019
Unfortunalely, even a million steps would not work :( I really appreciate your help though!
darova
darova on 14 Dec 2019
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
Cynthia Struijk on 14 Dec 2019
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 :)

Sign in to comment.

More Answers (0)

Sign in to answer this question.

Tags


Translated by