numerical solution of differential equations

[Ether]

Originally Posted by Hitchhiker 42

And the new sheets:

where did you get the revised accel formula?

[Hitchhiker 42]

Originally Posted by Ether

where did you get the revised accel formula?

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

[GeeTwo]

Originally Posted by Hitchhiker 42

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=Ae^3x/2

[Hitchhiker 42]

Originally Posted by GeeTwo

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

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

[Ether]

Originally Posted by Hitchhiker 42

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

In post14 you wrote:

Originally Posted by Hitchhiker 42

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]

Originally Posted by Ether

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]

Originally Posted by Hitchhiker 42

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]

Originally Posted by GeeTwo

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

Originally Posted by Hitchhiker 42

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

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

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]

Originally Posted by Hitchhiker 42

I've attached your spreadsheets modified for 1/4*x^3 + 1

Originally Posted by Ether

I am assuming what you actually meant by that is

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

Is that correct?

Originally Posted by Hitchhiker 42

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

Originally Posted by Hitchhiker 42

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]

Originally Posted by Ether

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]

Originally Posted by Hitchhiker 42

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. 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...rs/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]

Originally Posted by 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 + 2x³
Initial conditions: x(0) = 0, x'(0) = 1.

The analytic solution is simply:

x=tan t
x'=sec²t
x''=2 tant sec²t = 2tant + 2tan³t= 2x + 2x³

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, 10¹⁰⁰.

[Ether]

Originally Posted by Hitchhiker 42

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?