How to implement Forward Euler Method on a system of ODEs

Question

Cynthia Struijk 2019-12-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/496543-how-to-implement-forward-euler-method-on-a-system-of-odes

编辑： Cynthia Struijk 2019-12-14

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?

4 个评论
显示 2更早的评论隐藏 2更早的评论

Cynthia Struijk 2019-12-14

编辑：Cynthia Struijk 2019-12-14

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 2019-12-14

Where are the original equations?

Cynthia Struijk 2019-12-14

编辑：Cynthia Struijk 2019-12-14

 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 2019-12-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/496543-how-to-implement-forward-euler-method-on-a-system-of-odes#answer_406363

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 个评论
显示 3更早的评论隐藏 3更早的评论

darova 2019-12-14

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