I have discovered that in solving ordinary differential equations things can get very easy if we use micro-intervals (intervals divided by 10 000 or more) and very rough methods like Euler's , to get the same accuracy as when using complex methods.

ODEs of any order with initial conditions can be solved easely "using micro-intervals (intervals divide by more than a 1000) associated with Euler's method - definition of derivative -, can lead to great accuracy. I have a batch of scripts that can prove that. (study1, study2 ......). The simplicity is the great asset. Is this known ?

[Bjorn Gustavsson]

Yes. The error-characteristics of different numerical schemes for integration of ODEs are not unknown. See for example: Euler-errors

[Carlos Steinstrasser]

If I use micro-intervals plus Euler's method it's possible to obtain results as accurate as any other method with a dozen script lines. There are no visible truncation errors if you divide the interval by 10000 or even more. I have a master degree in Hydraulics and I just want to share this possibility with others because it opens a possibility to use the same procedure for partial differential equations. This is my line of research and my goal is to make things easier.

Thanks for your comment but maybe you did not get the point. All the math piruettes can be put aside if you use micro-intervals and the definition of derivatives.

[John D'Errico]

Is it known? Um, yeah. Any basic course in numerical methods will discuss things like this. At the same time, you need to consider that just using a tiny step size will often be really slow. As well, there are many issues with methods like Euler's method. You may find that even a tiny step size will not lead to convergence, especially on stiff systems. The problem is, there is no magic flag that pops up to alert you when Euler's method is giving you complete crapola because the system is stiff.

As methods go, Euler's method is a good one to teach to beginning students. Easy to follow. Easy to write code to implement. It is the equivalent of trapz. In fact, there is a very strong relationship between Euler's method and trapz. At the same time, neither is the right tool to use all the time. A good craftsman knows and understands the tools in their tool chest. But they also know when to not use a specific tool, and they would never use the same tool for all problems.

Learn the tools you have at your command, and learn when to use any given tool, and when not to do so too.

[Bjorn Gustavsson]

"...but maybe you did not get the point. All the math piruettes can be put aside if you use micro-intervals and the definition of derivatives..."

Yeah, I got the point all right. The point(s) you were supposed to read up on were how the error depends on the step-size (it is proportional to h^2 for your Euler) and its region of stability. One thing John points toward is that more efficient methods exist - that is methods that either give the same accuracy with larger step-size or better accuracy with the same step-size (or function-evaluations if you prefer that comparison), if I recall correctly the standard 4-th order R-K has an error-term proportional to h^5.

The math pirouettes are what you have mathematicians working on numerical analysis for! I am in no way any kind of expert on ODE/PDE-methods, but way back when we were taught that Euler was a really bad choise - poor stability, error-term only O(h^2) - better tools are available for as far as I understand every ODE-problem.

[Jim Riggs]

Carlos, you are absolutely right. You can use Euler's method rather than more sophisticated algorithms to solve your equations. It's only a mater of selecting the appropriate step size to control the error within your desired tolerence. The reason that this is true is because modern computers can represent numbers with a lot of digital precision. But this has not not always been the case.

I have a long career in the aerospace industry, virtually all of it developing and applying digital simulations as tools to solve design problems, so I am very familiar with the application of numerical integration tools. Runge-Kutta methods were oriinally developed as a means to solve equations while controlling the round-off error. In the early days, computers had much more limited digital precision available to represent floating point numbers, consequently, we spent a lot of time designing algorithms to control numerical error, and speed up execution by reducing the number of operations required. High-order numerical integration algorithms accomplish this by providing greater precision with fewer computations.

I am also a woodworker, and I know that any cut you can make using a table saw can also be done using a hand saw and a plane. Some would even argue that the quality of the resulting hand-cut is superior to that produced by the power tool. The main difference is the time and effort involved. Yes, there is a significant initial effort and expense in procuring and setting-up a table saw - adjusting the rip fence, miter gauge, and blade to be square and true. But once established, you can simply flip a switch and tear through cuts in almost no time with great accuracy. A long rip cut done with a hand saw, then trued with a hand plane might take an experienced craftsman 10 minutes, while the table saw can give a comprable result in about 30 seconds.

Anyone who owns a table saw will select it as the go-to tool over his hand saw whenever possible.

With modern computers and their large word length, Euler's method can give the same result as more complex, high-order algorithms, but will require several orders of magnitude more computation - requiring you to wait much longer for the answer, that's all. This can have a significant effect on one's productivity, so in the long run, a "professional" should equip himself with the most productive tools.

[John D'Errico]

Sorry, but a master's degree in hydraulics is not a degree in mathematics. It may just mean that you have not seen a sufficient variety of problems, so that for everything you have ever worked on, Euler's method would be fine. That is good for you, but that is not true in general. There are a huge variety of problems out there where your suggestion will fail to be a good choice, where Euler's method will not be stable, or where you would need such a small step size that it would take a huge amount of time to solve. Remember that some problems might involve complicated functions to evaluate, so that tens or hundreds of thouands of evals will be impossible to handle.

As I said above, Euler's method has no big flag that can tell you it is giving a poor solution.

Knowing what tool to use for different problems, as we have all said, is important. This is not just a question of mathematical pirouettes.

[Jim Riggs]

Here is a case study which I think is instructive. This illustrates the point that John is emphasizing:

Recently, I put together a simulation in order to compare the run time performance using Matlab compared to Python. The model is comprised simply of 2 bodies (Sun and earth) in mutual gravitayional attraction. With this simple model, the orbit of the earth about the sun has a very well defined theoretical behavior (i.e. it is a constant shape), therefore, any variation in the shape of the orbit is due to the numerical solver. I set up the model with the earth at it's average orbit radius (about 1.5e11 meters) and with an initial velocity computed to produce the requisite 365.25 day period. The momentum vector of the earth is computed, and the sun is initialized at the origin of the reference system, having an opposite momentum. (This establishes a fixed barycenter)

I run the model for a simulated time of 5 years. The solution error is computed from the drift in the orbit from it's initial starting point; i.e. how far the earth varies as it crosses the x-axis at the end of each orbit. The error is computed (in meters) as the max - min distance from the sun for the five orbits. Simulation runs were conducted for a matrix of time steps, and using Runge-Kutta integration methods (order 2 through 5). The Orbit dispersion error is shown in the following table:

So, for example, the entry labeled (1) in the table indicates that using a 4 hour time step and the 2nd order Runge-Kutta solver, the maximum dispersion in the earth orbit was 55,656.6 meters over the five year simulation period.

These errors show several things: First, the reduction in error with the reduction in the time step shows the classic "order of magnitude" improvement with each halving of the time step (for the RK2 and RK3 solutions). Note that as the time step goes from 4, 2, 1, .5, .2 the error goes (nominally) from 70,000, 7,000, 700, 70, 7 for RK2, and 30,000, 3,000, 300, 30, 3 for RK3. Second, there is a huge improvement going from RK3 to RK4, followed by no improvement at all going to RK5. This shows that the RK4 solver is a good fit to the model characteristics, and is the optimum choice for this problem.

Now, for comparison, I went back and added a forward Euler integration method to this model in order to compute the performance values for this method. Using a 4 hour time step, the orbit dispersion error was 2.08e10 meters (That's 13% of the orbit radius!). The following two charts show the reduction in the Euler method error as the time step is reduced (note that in the first chart , the time axis is in hours, where in the second chart the time axis is in seconds):

The lowest data point on the second chart had a time step of 3.6 seconds (0.001 hour) and shows an error of 6.7M meters, which is a little more than the radius of the earth. The execution time for this simulation run was 18 minutes. Using this data, noting that the error reduces linearly with the time step, I estimated what it would take to produce errors comprable to points (1), (2), and (3) in the first table, above. These are shown in red in the following table, next to the execution time for the Runge-Kutta methods.

Using the forward Euler method:

(1) In order to achieve an error of 55,656.6 meters the time step required is about 1/2 milisecond, and the estimated run time for the model is 37 hours. This compares to a run time of 0.56 seconds for the 2nd order Runge-Kutta implementation using a 4 hour time step.

(2) Reducing the error to 18,705.5 m would require a time step of .16 ms and an execution time of 108 hours, compared to 0.83 seconds using RK3 with a 4 hour time step.

(3) To achieve less than 7,000m error, a time step of 0.062 ms is required, and would run for about 292 hours, compared to 1.11 seconds using RK2 with a 2 hour time step.

This case study illustrates a case where Euler's method is a very bad option.

[Carlos Steinstrasser]

In a 64 bits computer, my experience with third order ODEs, using Euler method -takes a few seconds - after you have a smooth script. I'm talking about a simple script of about twenty lines, without any ready made functions. Try this, and you will see what I mean. Anyway you have done an excellent work.

[Steven Lord]

It is entirely possible that for the ODEs that you're solving as part of your work, using Euler's method is sufficient for your needs. Others may need to solve ODEs as part of their work for which Euler's method isn't sufficient. And there's nothing wrong with either of those statements.

MATLAB has eight different ordinary differential equations solvers whose names start with "ode" to support different accuracy / efficiency requirements. See the "Basic Solver Selection" table on this documentation page for a discussion of how the solvers differ and when you might want to use one over another.

That's not even counting the solvers available in Simulink, either.

[Carlos Steinstrasser]

In 1995 I bought "The Student Edition of Matlab" ,from that time on , I became a familiar with the Matlab developments in coding and all , meanwhile the personal computers have developed as everyone knows. I was amazed by the processing speed of the 64 bits PCs and had the idea of using micro-intervals plus very simple coding to achieve speed. The tips on coding by Matlab made me get a smooth processing and speed.

I'd like to share this idea with Matlab staff. You may use my "study1 ....to study8" scripts in my license. It's always the same structure with different ODEs.

The back and forth pirouettes of the past (Runge-Kutta 4 etc) maybe a problem for speed in the knew Pcs.

Thanks.

[James Tursa]

+1 @ Jim Riggs for good example ... Euler Method will systematically always overshoot and build up error. And even with very tiny step sizes, I think one would have a hard time getting the accuracy with Euler Method because floating point word size effects would start to dominate.

[Carlos Steinstrasser]

This is not a guess. See my scripts and try by yourself. I don't use ready-made Euler method.

[Carlos Steinstrasser]

%y''=xy'-4y   y(0)=3.0 e y'(0)=0.0  exact solution(y=x^4-6x^2+3)
%y'=z
%y''= z'    
%z'=xz-4y

clear
dx=0.00001;
ilim=50000;  %dx*ilim=0.5 constante xlim
x=0.0:dx:0.5;

y(1)=3.0;
z(1)=0.0;
x(1)=0.0;

for i=1:ilim
    y(i+1)=y(i)+dx*z(i);
    z(i+1)=z(i)+dx*((z(i)*x(i))-4.0*y(i));
    
end
plot(x,y)
hold on
xx=0.0:0.05:0.5;
for ii=1:11
    yes(ii)=(xx(ii)^4)-(6.0*xx(ii)^2)+3.0;
end

plot(xx,yes,'r.')
hold off

%yes(1:ilim+1)=(x(1:ilim+1).^4)-6.0*(x(1:ilim+1).^2)+3.0;

[James Tursa]

@Carlos: From reading the various comments, it does not appear to me that anyone is saying Euler's Method with tiny step sizes does not work well for the particular problems you have been investigating. Maybe it does. What we are saying is that Euler's Method is known to be bad for lots of other problems (see Jim Riggs example for one, or stiff problems mentioned by John D'Errico), and can lead to horrendous run times even in cases where it might achieve the accuracy desired. So as a general technique as you seem to be proposing, I don't think you will find anyone in this forum in agreement with you. And, as others have already mentioned, this idea is nothing new and studies in step sizes vs methods have been around for a long time.

[Carlos Steinstrasser]

Give me a stiff problem with clear initial conditions and let me try to solve it using my naive ideas. The ODEs or system shall be clearly stated and also the solution I'm supposed to find.

[Steven Lord]

Try the example in the Numerical Stability section of the Wikipedia page for Euler's method.

As a second example, try the flame propagation example in this post on Cleve's blog.

[Carlos Steinstrasser]

Thanks for the example.

Please type the 9 lines of the script and see what I mean.

%y'=-2.3*y   y(0)=1.0   exact solution(y=exp(-2.3x)
clear
dx=0.001;
ilim=5000;  %dx*ilim=5.0 constante xlim
x=0.0:dx:5.0;
y(1)=1.0;
for i=1:ilim
    y(i+1)=y(i)-2.3*y(i)*dx;
end
plot(x,y)

%The instability is only with dx=1.0; micro-intervals and Euler is the answer

%run this script y(5001)=>y(5)=9.9968*10^(-6)

[Bjorn Gustavsson]

In my field of work I often need to solve the equations of motion for charged particles in E and B-fields. The simplest test case is this ODE:

% dx/dt = vx;
% dy/dt = vy;
% dvx/dt = vy;
% dvy/dt = -vx;
odevxB = @(t,rv) [rv(3);rv(4);rv(4);-rv(3)];

That is the 2-D equation of motion for a charged particle in a homogenous magnetic field, in normalized units. Integrate that one from 0 to , with initial condition:

rv0 = [-1,0,0,1];

For that the first numerical test is to verify that the ODE-integration-scheme conserves the constants of motion, in this case the kinetic energy of the particle should be conserved (the force perpendicular to the velocity) and therefore also the gyro-radius. So when done, check these plots:

subplot(211),
% with x-component in the first and y-component in second 
% column this should be 1 (gyro-radius squared)
plot((rv(:,1)-mean(rv(:,1))).^2 + (rv(:,2)-mean(rv(:,2))).^2)
subplot(2,1,2)
% velocities in 3rd and 4th columns, this should be proportional 
% to the kinetic energy
plot(rv(:,3).^2 + rv(:,4).^2)

Does your explicit Euler-integration conserve what it should?

[Carlos Steinstrasser]

@Steven Lord: This is about the flame problem:

%y'=(y^2.0)-(y^3)   y(0)=0.02   
%y'=(y^2.0)*(1-y)
clear
dx=0.002;
ilim=50000;  %dx*ilim=5.0 constante xlim
x=0.0:dx:100.0;
y(1)=0.02;
for i=1:ilim
    y(i+1)=y(i)+dx*(y(i)^2.0)*(1.0-y(i));
end
plot(x,y)

These few lines solve the problem. I just wanted to show that ODE can be solved with simple methods. Thanks to all the people that wasted their time following my ideas.

Matlab has been a great help in my academic life.

[Jim Riggs]

Carlos,

I coded your equation, above, using Euler method and 100,000 calculation steps:

I also computed the solution using the Runge-Kutta methods, order 2 through 5.

The exact solution 0 <= x <= 0.5 is 1.5625

For the Euler method, using 100,000 calculation steps, the solution value is 1.56251011, resulting in an error of 0.00001011. (this ammounts to 6.47e-04 % error) I then adjusted the number of steps required to obtain this same accuracy (or better) using the Runge-Kutta algorithms. The table shows the results.

As you can see that there is a vast difference in efficiency of the Runge-Kutta methods compared to Euler. This is the reason that they are so popular. For high-order equations, using a high-order solver significantly reduces the amount of computation required to obtain a comprable result. This means faster performance. For this example, the RK2 solver is 125x faster than Euler, and RK5 is nearly 1000x faster.

Below is a sample of my code. It is not that much more work to utilize all the power available in multi-pass solvers. The most convenient way to apply this is to define a "state vector" which contains all of the states to be integrated, and the "state derivative vector", which is the vector of the derivatives for each state. The solver function "RK_integration" is coded to operate on vector inputs.

clear xsave
% select solver options
RKorder = 4;
outputrate = 1;  % =1 saves every point,  =2 saves every other point, etc.
nsteps = 100000;    % number of calculation steps

% Initialize the user function values and state vector
y = 3.0;
z = 0.0;
x_start = 0.0;
xmax = 0.5;

S(1) = y;  % state vector
S(2) = z;

% initialize remaining parameers
cnt=1;      % count the number of integration steps
nout=1;     % count the number of output points
firstpass = true;
x = x_start;
dx = (xmax-x_start)/nsteps;   % calculate the step size

% add the initial values to xsave 
xsave(nout,1) = x;
xsave(nout,2:3)=S;

%============ Start Solution loop ============
while x <= xmax
    DS(1) = S(2);  % DS is the state derivative vector, S is the state vector
    DS(2) = x * S(2) - 4*S(1);
% note that the integration function updates the state vector, S, and the value of x
    [S,x,trigger] = RK_integration(S,DS,dx,x,RKorder,firstpass);
    firstpass = false;
    
% log the output at the end of each integration step (trigger = true)
    if trigger
        cnt=cnt+1;   % count the number of outputs
        if mod(cnt,outputrate)==0  % this allows me to skip updates
            nout=nout+1;         % count number of saved values
            xsave(nout,1)=x;     % save x
            xsave(nout,2:3)=S;   % save the state vector
        end
    end
end

With this structure, the integration function controls the value of x, and performs as many itterations as required by the selected algorithm. The logical "trigger" allows the integration function to indicate when it has completed enough evaluations and is ready to move on to the next calculation interval. For example, setting RKorder=4 requires the integration algorithm to perform 4 evaluations of the function, using updated values of x and the state vector S along the way. (for this exercise, I modified the RK_integration function to add the Euler integration option for RKorder=1).

Seeing how easy this is to implement, why would anybody confine themselves to the Euler method?

For my 2-body orbit model, reference in my earlier comment, simply replace the two lines of code calculating DS with a function call to compute the state derivatives. The n-body subroutine is only 57 lines of code, and can accomodate any number of bodies. There are 6 states for each body (3 positions and 3 velocities), so the 2-body state vector has 12 states. In fact, you can replace the two lines that calculate DS with an entire simulation, consisting of thousands of lines of code.

The only other change required in the above script is to initialize the state vector per the requirements of the particular problem (and adjust the size of the xsave array)

It's really easy to set up the solver for a new user function, so why not use all the power of high-order integration methods? They are truly far superior, and much faster.

[Carlos Steinstrasser]

Maybe you are right. I'd just like to raise a flag that what was a great idea in the past may have to be reconsidered because of the new computers. In a very basic work, and I'm not an specialist in anything, I realised that every ODE I have tried so far can be solved with micro intervals and Euler's method . Of course they can be solved with Runge Kutta methods or other more recent ones. I thought this could be useful for new students because of the simplicity of the ten script lines. I'm thankful for all the attention I had from you all.This is all for reflexion and to face the new chalenges that certainly will show up.

Thanks Jim Riggs.

[Bjorn Gustavsson]

Carlos, did you try the equations of motion I suggested above?

[Carlos Steinstrasser]

@Bjorn: No I didn't, because I need to understand completely the model and what I'm looking for. Maxwell equations are not "my thing". Otherwise, the system of first order equations are very basic in appearance and Euler certainly can be used. See the Jim Riggs's explanation above.

[Bruno Luong]

"I'd just like to raise a flag that what was a great idea in the past may have to be reconsidered because of the new computers. "

New computer? The double precision floating point 64 bits exists almost at the born of the computers, at least already in 1953 with FORTRAN, more than 50 years ago.

Not sure about the history, but it's not a surprise me that success of the moon landing in 1969 is greatly due to RK methods.

And the study of errors of RK methods mainly focus on the high order approximation of the solutions and little related to the truncation AFAIK.

[Bjorn Gustavsson]

OK, I can help explain the equation. As you've understood it is a simple system of first order equations. It is, however, not at al related to Maxwell's equations, it is the Newtonian equation of motion:

Further, you are correct in that an explicit Euler-method can be used. That is not the question here. The question is whether Euler-methods are suitable for this problem. For these kind of clean physical problems we know that there are constants of motion that should be conserved. Here the total energy is fixed - in this case only the kinetic energy of the single particle, i.e. should be constant, the particle should neither gain nor loose energy, mening the speed should be constant. This requirement is something that is in no way a given for general-purpose ODE-methods, Euler, R-K or other schemes, but a complementary abd completely different physcis-based requirement. For this ODE-schemes like the Størmer-Verlet or Boris-mover are preferable. So this question falls outside from Jim Rigg's eplanation above.

The equation as implemented above are in normalized units such that the gyro-period is , the first two components are for the time-derivative of particle position (x and y) and the two last components are the corresponding accellerations.

Answer 1

[Jan]

Yes, this is known for a long time already. As Euler has known already and as you learn in the first year of numerical maths, the local truncation error of the Euler method is proportional to h^2. This means that reducing the step size will decrease the local truncation error also. In addition this will increase the number of function evaluations and the accumulated rounding errors. For many real world problems the later effect will dominate the run time and the total error of the computed trajectory. E.g. the ODEs to simulate the deformation of a car during a crash can include thousands or millions of variables and a function evaluation can take seconds or minutes. Then it is a bad idea to reduce the step size by a factor of 10'000. In addition this can let the accumulated rounding error explode. In consequence the Euler method is not used for scientific applications, but higher order methods with a smart step size control to find the minimal total error.