**You are now following this question**

- You will see updates in your activity feed.
- You may receive emails, depending on your notification preferences.

43 views (last 30 days)

Show older comments

Bjorn Gustavsson
on 15 Jul 2019

Carlos Steinstrasser
on 16 Jul 2019

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
on 16 Jul 2019

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
on 16 Jul 2019

"...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
on 16 Jul 2019

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
on 16 Jul 2019

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
on 17 Jul 2019

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
on 17 Jul 2019

Steven Lord
on 17 Jul 2019

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.

Carlos Steinstrasser
on 17 Jul 2019

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
on 17 Jul 2019

Carlos Steinstrasser
on 17 Jul 2019

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

Carlos Steinstrasser
on 17 Jul 2019

%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
on 17 Jul 2019

Carlos Steinstrasser
on 18 Jul 2019

Carlos Steinstrasser
on 18 Jul 2019

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
on 18 Jul 2019

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
on 18 Jul 2019

@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
on 18 Jul 2019

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
on 19 Jul 2019

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.

Carlos Steinstrasser
on 19 Jul 2019

Bruno Luong
on 19 Jul 2019

"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
on 19 Jul 2019

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.

Jan
on 19 Jul 2019

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

Select web siteYou can also select a web site from the following list:

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

- América Latina (Español)
- Canada (English)
- United States (English)

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)