How to implement Forward Euler Method on a system of ODEs

2019 12 13
1 个回答
回答已采纳
更新时间：2019 12 14
14 次查看（30 天）

0 个投票

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更早的评论

darova
darova 2019-12-14
Yes, it's possible. Do you have any attempts?
Cynthia Struijk
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
darova 2019-12-14
Where are the original equations?
Cynthia Struijk
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] ];

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

darova
darova 2019-12-14

0 个投票

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更早的评论

Cynthia Struijk
Cynthia Struijk 2019-12-14
编辑：Cynthia Struijk 2019-12-14
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
darova 2019-12-14
Try to increase number of steps
n = 1000; % number of steps
Cynthia Struijk
Cynthia Struijk 2019-12-14
Unfortunalely, even a million steps would not work :( I really appreciate your help though!
darova
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
Cynthia Struijk
Cynthia Struijk 2019-12-14
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 :)

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by