numerical solution of differential equations

[Ether]

consider the function x = cos(t)

the first derivative is x' = -sin(t)

and the second derivative is x" = -cos(t)

so we have the differential equation x" = -x

and x = cos(t) is the analytical (true) solution to that differential equation
with initial conditions x=1 and x'=0 at t=0.


Now turn things around.


suppose we are given the differential equation x" = -x',

with initial conditions x=1 and x'=0 at t=0,

and we want to plot x and x' vs t,

but we don't know how to find the analytical solution,

so we decide to numerically integrate it using the Euler method:

x'[n+1] = x'[n] + x"[n]*dt

x[n+1] = x[n] + x'[n]*dt

x"[n+1] = -x[n+1]

See attached spreadsheet Euler.XLS to see what happens. Yikes.

Columns Xa and Va are the analytical (true) solutions for position and velocity.

Columns Xe and Ve are the Euler Method numerical solutions for position and velocity.



Now instead of using the Euler method, use the Midpoint method:

Vmid = x'[n] + x"[n]*dt/2

Xmid = x[n] + Vmid*dt/2

Amid = -Xmid

x[n+1] = x[n] + Vmid*dt

x'[n+1] = x'[n] + Amid*dt

x"[n+1] = -x[n+1]


See spreadsheet Midpoint.XLS

Columns Xa and Va are the analytical (true) solutions for position and velocity.

Columns Xm and Vm are the Midpoint Method numerical solutions for position and velocity.

Columns Vmid, Xmid, and Amid are the extra columns needed for the Midpoint Method.

Notice that even though the Midpoint method requires 3 additional columns,
you can double the step size dt,
so the computation is just as fast as Euler,
but with far better accuracy.

[Richard Wallace]

Hi Russ!

I see you are bored, or maybe just antsy for the 2017 FIRST Robotics Competition kickoff. Either way, this is an interesting topic in numerical methods.

Question for an interested student (someone much younger than me):

Should we expect, in general, that Euler integration will be well suited to approximate monotonic systems, while Midpoint integration gives better results for periodic systems? If so, why?

[Ether]

maybe just antsy for the 2017 FIRST Robotics Competition kickoff.

Yes I am very much looking forward to the Kickoff.

And I have some long-overdue hardware to return to you :)

[Hitchhiker 42]

Hi Russ!
Should we expect, in general, that Euler integration will be well suited to approximate monotonic systems, while Midpoint integration gives better results for periodic systems? If so, why?

As I understand it, the main difference between the Euler method and the Midpoint method is that the midpoint method takes the slope by connecting a point behind and a point ahead of the given point. The Euler method just takes the slope at that point at extrapolates for that step. Please, do correct me if I'm wrong.

So, a monotonic function only increases (or decreases). It seems Euler integration would be better suited (read: more accurate) because the midpoint method depends on points behind the given point that's being calculated, which will tend to keep the slope smaller than it should be, whereas because the function doesn't tend to change direction (up or down) as much (it can only go one direction - monotonic), approximating ahead will tend to be better than approximating while taking into account behind the current point as well.

[GeeTwo]

What surprises me more than the gain in amplitude for Euler (which is pretty easy to guess if you consider what happens to energy at different points) is the excellent prediction of the period. I'll have to give this a look.

I was able to do a version without the two extra columns that tracked pretty closely, using the parabolic formula for constant acceleration to calculate the next position, and the average acceleration assuming constant jerk (x''') to calculate the next velocity.

x[n+1] = x[n] + dt*(x'[n] + x''[n]*dt/2)
x''[n+1] = -x[n+1]
x'[n+1] = x'[n] + dt*(x''[n] + x''[n+1])/2

[Ether]

As I understand it,...the midpoint method takes the slope by connecting a point behind and a point ahead of the given point...Please, do correct me if I'm wrong.

Where did you come by this understanding?

It seems Euler integration would be better suited [for monotonic function]...

You could easily modify the XLS I posted to test this hypothesis :)

[Ether]

Unless I was careless with the algebra (it happens when I'm tired),
Steve's catapult (https://www.chiefdelphi.com/forums/showthread.php?t=152739) can be modeled with an ODE of the form

θ'' = k1 + k2∙cos(θ) + k3∙θ'

Attached is an Octave script that uses Octave's built-in ODE solver "lsode"
to numerically integrate arbitrary ODEs of the form x'' = f(t,x,x')

[sspoldi]

Hillbilly solution:
one forward Euler integration, one backward Euler integration, less typing and good enough for government work.

Surprising how often that works...

Cheers,
Steve.

P.S. The equation (for the catapult) should be something like θ" = K1∙(K2 - θ'), θ is just along for the ride.

[GeeTwo]

Hillbilly solution:
one forward Euler integration, one backward Euler integration, less typing and good enough for government work.

Surprising how often that works...

That is:

x'[n+1] = x'[n] + dt*x''[n] //backward looking
x[n+1] = x[n] + dt*x'[n+1] //forward looking
x''[n+1] = -x[n+1]

It looks OK on amplitude, but overpredicted the resonant frequency by .. almost half a part per thousand. Certainly good enough for FRC.

P.S. The equation (for the catapult) should be something like θ" = K1∙(K2 - θ'), θ is just along for the ride.
The cosθ term is gravity acting on the boulder (and lever arm).

[sspoldi]

The cosθ term is gravity acting on the boulder (and lever arm).

Yea, I'd like to say we just ignore stuff like that, but sometimes it's a real factor. In 2014 we had a hammer with a 3 pound head on a 1 foot arm, gravity definitely made a difference.

Since we typically don't have a lot of time (who does), I like to get the kids to do a simple model up front, and then we do some system id and fit the actual robot behavior to a model. This way we can tune control systems quickly, and it gives them a chance to do some data based optimization in addition to a little physics up front.

Cheers,
Steve.

[Ether]

one forward Euler integration, one backward Euler integration

Vn+1 = Vn + An∙dt;
Xn+1 = Xn + Vn+1∙dt;

Maybe provides some insight: the above is algebraically equivalent to

Vn+1 = Vn + An∙dt;
Xn+1 = Xn + dt∙(Vn+Vn+1)/2 + ½∙An∙dt2

[GeeTwo]

Maybe provides some insight: the above is algebraically equivalent to

Vn+1 = Vn + An∙dt;
Xn+1 = Xn + dt∙(Vn+Vn+1)/2 + ½∙An∙dt2

..so counting the constant acceleration term twice in the position calculation mostly offsets not counting the jerk in the velocity calculation..at least in this case.

[Ether]

θ'' = k1 + k2∙cos(θ) + k3∙θ'

derivation of k1 k2 k3

[Hitchhiker 42]

You could easily modify the XLS I posted to test this hypothesis :)

I've attached your spreadsheets modified for 1/4*x^3 + 1.
It seems that the Euler method does approximate it better.

[Ether]

@Mark: Where did you get the accel formula for columns D and G.

[Hitchhiker 42]

@Mark: Where did you get the accel formula for columns D and G.




Not sure. Must have been really tired :) . I did that over, and also added an error column to show the difference between Xa and the predicted X (Xm or Xe).
I set dt = 0.01 s and compared the errors at t = 1, 2, 4 s

For 1s:
Midpt error: 0.598
Euler error: 0.585

For 2s:
Midpt error: 2.834
Euler error: 2.748

For 4s:
Midpt error: 50.070
Euler error: 48.113

Seems that, at least for this function, the Euler method holds up better. Tomorrow, I'll try it with a function where the second derivative isn't always positive on the interval I'm testing.

[Hitchhiker 42]

And the new sheets:

[Ether]

And the new sheets:

where did you get the revised accel formula?

[Hitchhiker 42]

where did you get the revised accel formula?





I took the derivative of the actual velocity function. Is this not correct?

[GeeTwo]

I took the derivative of the actual velocity function. Is this not correct?

No, as I read it you're setting x''=3x/2. This would be solved by a function of the form x=Ae3x/2

[Hitchhiker 42]

No, as I read it you're setting x''=3x/2. This would be solved by a function of the form x=Ae3x/2

Could you explain what this means? I'm not sure I understand...

[Ether]

Could you explain what this means? I'm not sure I understand...

In post14 you wrote:
I've attached your spreadsheets modified for 1/4*x^3 + 1.

I am assuming what you actually meant by that is

x(t) = 1/4*t^3 + 1

Is that correct?

If so, then in order to investigate the performance of Euler vs Midpoint numerical solution (in the context of this thread) you need to convert that to an initial value problem (IVP) consisting of a differential equation (involving only x, x', and x'') plus initial values for x(t) and x'(t) at some point t=to.

[Hitchhiker 42]

In post14 you wrote:


I am assuming what you actually meant by that is

x(t) = 1/4*t^3 + 1

Is that correct?




Yes, sorry about that. That was what I meant.

[Ether]

Yes, sorry about that. That was what I meant.

OK.

x(t) = 1/4*t^3 + 1 is the analytical solution to the Initial Value Problem

a(t) = sqrt(3*v(t))

with initial values*

x(0.01) = 1.00000025

and

v(0.01) = 7.5e-5

As you can see from the spreadsheet, the Midpoint method gives much better results for that IVP than Forward Euler and Forward/Backward Euler.


* selected to avoid the stationary point at t=0

[GeeTwo]

No, as I read it you're setting x''=3x/2. This would be solved by a function of the form x=Ae3x/2

Could you explain what this means? I'm not sure I understand...

Sorry, that should have been x''=Ae√3t/√2:o

if x=Ae√3t/√2,
then x'=√3Ae√3t/√2/√2 (chain rule)
and x''=3Ae√3t/√2/2 (chain rule again) = 3x/2.

Edit: I knew there was more to it, and just figured out the other part:

x''=3x/2 ==> x= Ae√3t/√2 + Be-√3t/√2

[Ether]

I've attached your spreadsheets modified for 1/4*x^3 + 1
I am assuming what you actually meant by that is

x(t) = 1/4*t^3 + 1

Is that correct?
Yes, sorry about that. That was what I meant.

I took the derivative of the actual velocity function. Is this not correct?

Yes, but you got x and t confused.

given:

x(t) = 1/4*t^3 + 1

... take time derivative of x(t) to get:

x'(t) = (3/4)*t^2

...take time derivative of x'(t) to get:

x''(t) = (3/2)*t


Now solve for the differential equation:

solve x'(t) for t:

x'(t) = (3/4)*t^2 => t=sqrt(4x'(t)/3)

... and then substitute for t in x''(t):

x''(t) = (3/2)*t = (3/2)*sqrt(4x'(t)/3) = sqrt(3*x'(t))

[Hitchhiker 42]

Yes, but you got x and t confused.

given:

x(t) = 1/4*t^3 + 1

... take time derivative of x(t) to get:

x'(t) = (3/4)*t^2

...take time derivative of x'(t) to get:

x''(t) = (3/2)*t


Now solve for the differential equation:

solve x'(t) for t:

x'(t) = (3/4)*t^2 => t=sqrt(4x'(t)/3)

... and then substitute for t in x''(t):

x''(t) = (3/2)*t = (3/2)*sqrt(4x'(t)/3) = sqrt(3*x'(t))




I see... so that column is in terms of x', not t. I've updated the spreadsheet to reflect that and it seems to show that the Euler method does better.
dt = 0.01s, error measured at t=3s.

Euler error: 0.2061
Midpoint error: 5.0828

This, I presume, is for the reason I mentioned in my first post.

[Ether]

I see... so that column is in terms of x', not t. I've updated the spreadsheet to reflect that and it seems to show that the Euler method does better.

You're still struggling with this, but I can't point out your error since you didn't post your spreadsheet.

I think you may have overlooked my earlier post #24 (https://www.chiefdelphi.com/forums/showpost.php?p=1623411&postcount=24). The attachment shows how to set up the formulas.

The Midpoint method is clearly better than either Forward/Forward Euler or Forward/Backward Euler for this IVP.

[Hitchhiker 42]

I redid the midpoint spreadsheet - now at t=9.96s with dt=0.01s, the error is only 0.96. I replaced the Am column with the Am = sqrt(3*x').

Now to search for a possible explanation...
I did find this http://www.physics.drexel.edu/~steve/Courses/Comp_Phys/Integrators/simple.html

My guess would be that for this function, which, on the interval being examined, is concave up always (x'' > 0), the Euler method uses a point from farther back than the midpoint method to determine the slope (beginning of the interval slope vs. mid-interval slope). As it's concave up, the Euler method lags behind the midpoint method in growth, and so grows slower than the actual function by more than the midpoint method.

Now to test a method where it uses the end of the interval's slope...

[Hitchhiker 42]

Now to test a method where it uses the end of the interval's slope...

Drat... I tried it using the slope at the end of the interval instead of at the beginning, and it turned out to have a tad bit more error compared to the actual function.

[GeeTwo]

The key to efficiently numerically solving this sort of differential equation is to understand the size/scale of the various derivatives. If the second derivative is consistently small, a linear (simple Euler) integration can provide a decent prediction at a reasonably large range step. As the second derivative increases, the range step must be reduced and/or the second derivative must be included in the calculation for each range step; midpoint is one way to do this. If the third derivative is large, you need to make the steps even smaller! Fortunately, the importance of the nth derivative is scaled by 1/n!, that is, the reciprocal of n factorial (google Taylor Expansion if you aren't familiar with it).

If you want an extreme case to work with that has a simple analytic solution, try this one:

x''(t) = 2x + 2x3
Initial conditions: x(0) = 0, x'(0) = 1.

The analytic solution is simply:
x=tan t
x'=sec2t
x''=2 tant sec2t = 2tant + 2tan3t= 2x + 2x3

Using a fixed step size in t, and a finite number of terms in the Taylor expansion with Euler-like integration, the numeric solution will never reach infinity, whereas the analytic solution reaches infinity at pi/2. Try to find solutions and step sizes which suitably predict pi/2, based on x being greater than, say, 10100.

[Ether]

I tried it using the slope at the end of the interval instead of at the beginning.

It's not clear what the point of that spreadsheet is. You are using the analytical solution for x' to compute an Euler step.

You weren't intending to compare that result to the Midpoint method, were you?