Using Euler's method as ODE-solver in Matlab

65 views (last 30 days)
Katara So
Katara So on 3 Mar 2021
Edited: Jan on 4 Mar 2021
I want to implement an ODE-solver to solve a set of given differential equations to give the same figure as the one below and also study the accuracy order of the method.
My differential equations are:
f = @(t,y) [0.272 - 0.00136*y(1) - 0.00027*y(1)*y(4);
2.7e-05*y(1)*y(4) - 1.36e-3*y(2) - 3.6e-2*y(2);
0.243e-3*y(1)*y(4) + 3.6e-2*y(2) - 0.33*y(3);
100*y(3) - 2*y(4)];
where
t0=0; tend=120; y0=[200;0;0;4*10^-7]
How can I use Euler's method?
I found an example using the same differential eqs. as the ones I was given (which are from an article) I will post it below. However, I don't quite understand what is being done in the code. How was the step size determined to be h=1/18? What does Euler's method in the while sling actually do? And from where did the f1 function come from, what does it calculate? Furthur, how can the accuracy order be determined from this? I have never done anything like this before so all help is truly appreciated!
% y = [R;L;E;V] = [y(1);y2;y(3);y(4)]
y0 = [200;0;0;4*10^-7]
f = @(t,y) [0.272 - 0.00136*y(1) - 0.00027*y(1)*y(4);
2.7e-05*y(1)*y(4) - 1.36e-3*y(2) - 3.6e-2*y(2);
0.243e-3*y(1)*y(4) + 3.6e-2*y(2) - 0.33*y(3);
100*y(3) - 2*y(4)];
t = 0;
ts = 120;
h = 1/18;
y = y0;
yout = y0;
while t < ts
y = y + h*f(t,y);
t = t + h;
yout = [yout y];
end
v = (0:1/18:120)
f1 = 800 + yout(1,:) + yout(2,:) + yout(3,:)

Answers (1)

Jan
Jan on 3 Mar 2021
How was the step size determined to be h=1/18?
The author of the code wanted it to be 1/18.
What does Euler's method in the while sling actually do?
It computes the derivative at the current point by evaluating f() and approximates the curve by going a step in this direction as a straight line.
And from where did the f1 function come from, what does it calculate?
If you mean f(), this function was chosen by the author of the code. All we know as readers is the formula. We cannot have the faintest idea if this is the price of bananas or the trajectory of a ball on a curved surface.
Furthur, how can the accuracy order be determined from this?
A simple investigation is to reduce the stepsize h and compare the resulting trajectory. In numerical maths you learn a lot about the accuracy of the Euler method.
  2 Comments
Katara So
Katara So on 3 Mar 2021
Thank you!
If I want to plot it, changing the stepsize h gives me an error that v and f1 are not of same length for some h:s while it works for others. For ex, h =1/18 works, but h=1/20 does not and h=1/25 works. Why is that? How did the author know that h = 1/18 would work?
Jan
Jan on 4 Mar 2021
Edited: Jan on 4 Mar 2021
"gives me an error that v and f1" - Please do not rephrase the error message, but post a copy of the complete message. Then the readers have better chances to understand, what's going on.
The solution is easy:
t = 0;
ts = 120;
h = 1/18;
x = [];
while t < ts
t = t + h;
x = [x, 1];
end
v = (0:1/18:120)
size(x)
size(v)
Note the "t < ts": If ts is a multiple of h, the last step reaches ts exactly (except for rounding errors). With "t <= ts" this would create an x with the same number of elements as 0:h:120. But with "t < ts" the element for t=120 is omitted in the WHILE loop.
If 120 is not a multiple of h, this effect does not appear, see:
1:4:10
For this reason the WHILE loop is a bad choice for the Euler method, because the final time is not reached under some conditions. Fixing the time span and h cannot reach the final time point. Either append a final step after the while loop to reach ts or define the number of steps instead of h. Then:
t = linspace(t0, ts, NumberOfSteps)
h = t(2) - t(1);

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!