numerical solution of differential equations
 

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.

Re: numerical solution of differential equations
 

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?

Re: numerical solution of differential equations
 

Yes I am very much looking forward to the Kickoff.

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

Re: numerical solution of differential equations
 

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.

Re: numerical solution of differential equations
 

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.

Re: numerical solution of differential equations
 

Where did you come by this understanding?

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

Re: numerical solution of differential equations
 

Unless I was careless with the algebra (it happens when I'm tired),
Steve's catapult can be modeled with an ODE of the form

θ'' = k₁ + k₂∙cos(θ) + k₃∙θ'

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')

Re: numerical solution of differential equations
 

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.

Re: numerical solution of differential equations
 

That is:

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

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

Re: numerical solution of differential equations
 

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.

Re: numerical solution of differential equations
 

V_n+1 = V_n + A_n∙dt;
X_n+1 = X_n + V_n+1∙dt;

Maybe provides some insight: the above is algebraically equivalent to

V_n+1 = V_n + A_n∙dt;
X_n+1 = X_n + dt∙(V_n+V_n+1)/2 + ½∙A_n∙dt²

Re: numerical solution of differential equations
 

..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.

Re: numerical solution of differential equations
 

derivation of k₁ k₂ k₃

Re: numerical solution of differential equations
 

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

Re: numerical solution of differential equations
 

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