New algorithms for numerical integration and ODE solutions

[Namir]

Hi All,

I have posted on my web page two new algorithms for numerical integration and ODE solutions. Both methods share the same strategy that allow you to control the calculations at approximately a regular pace for the integration with the dependent variable Y instead of the typical approach of controlling the integration pace with the independent variable X. I have included Excel files that show how the algorithms work using VBA code.

Enjoy!

Namir

[Dan W]

Hi Namir

Your SONI algorithm uses this:
h = 2 * |YMaxDelta / f’(X1)|

When I used the function EXP(-X)*(SIN(X))^2, this function evaluates to zero at x=0, and also the derivative is zero at x=0. This should cause the equation for h to blow up.

It doesn't in your example because your MyDx function is an approximation. However if an improved MyDx function is used that generates a number much closer to zero, h ends up being very large.

In this case the next DO loop evaluates the function at a large X2, which coincidentally is also close to zero causing it to exit the DO.

I think some protection on h near zero slope and when the slope is rapidly changing would be a valuable improvement.

[Namir]

Dan,

I have coded the approximation for the derivative to be:

f'(x) = (f(x+h) - f(x)) / h

And not the more accurate form:

f'(x) = (f(x+h) - f(x - h)) / 2h

To avoid the case when f(x+h)=f(x-h) around x=0, which would lead to an error-geneating zero slope.

Namir

[Namir]

I went ahead and updated the Excel VBA code and the documents to reflect handling zero and small slopes.

Namir

[Dan W]

When I ran your mod using Gauss2 on Sheet10, I got a result that was not very accurate. the result was something like 0.69.

Because of the small slope at zero, the first value of h was 1000, then Y2=0 and the calculation of h is exited. X2 is then reduced to 5 (XMax), then the Gauss algorithm calculates the integral in one interval from 0 to 5.

I recommend modifying the calculation of h as follows:

XMaxDelta = (Xmax - X1) / 5
' Put an upper limit on the interval size

h = 2 * Abs(YMaxDelta / MyDx(sFx, X1))
If (h > XMaxDelta) Then h = XMaxDelta
For i = 1 To 100 ' prevent infinite loop (may not be necessary)

h = h / 2
X2 = X1 + h
Y2 = MyFx(sFx, X2)
If (Abs(Y2 - Y1) <= YMaxDelta) Then Exit For

Next i

A better way might be to remove EPSILON from MyDx (allow that function to be as accurate as possible) and calculate 1/h first to avoid a divide by zero, then limit 1/h.

[Namir]

Thank you Dan for your feedback. I will add an Xmax value. The ODE version uses steps for X to display intermediate Y values. These steps act as a type of Xmax value.

Namir

[Namir]

I updated the Slope-Oriented Numerical algorithms on my web site.

Namir

[Wes Loewer]

(03-13-2014 10:44 AM)Namir Wrote:  I have posted on my web page two new algorithms for numerical integration

Namir,

Thank you for the articles.

I have a question about your Romberg-Simpson Method. (I've seen it mentioned in passing that this could be done, but I've never actually seen it implemented.) Your algorithm uses:

R(n,m) = (1/4^m-1)*(4^m*R(n,m-1)-R(n-1,m-1))

for both the Trapezoid and Simpson methods.

The derivation that I've seen for this recursive formula was based on the fact that the Trapezoid Error is proportional to 1/n^2. I'm not certain about this, but since the Simpson Error is proportional to 1/n^4, it seems like the recursive formula for Simpson's would be different. I could be totally wrong, but using a similar derivation, I'm coming up with:

R(n,m) = (1/16^m-1)*(16^m*R(n,m-1)-R(n-1,m-1))

I tried a couple of examples on a spreadsheet and this second formula seems to converge faster. Try it in your VB code and compare the two.

Please correct me if I'm off on this.

-Wes L

[Namir]

?????????????

In the Slope-Oriented intergration aritcles I do not use the Romberg method as one of the basic integration methods called by my method.

Are you looking at another article on my web page?

Namir

[Wes Loewer]

(03-26-2014 09:26 PM)Namir Wrote:  Are you looking at another article on my web page?

Yes, one of your older articles, "A New Face of Romberg Integration," caught my eye and I didn't realize till after I posted that this wasn't one of the "new algorithms" you mentioned.

-Wes L

[Namir]

Wes,

I try to use 16 instead of 4, but it seems I failed miserably in my updates. One function gave more error while another one cause run-tine overflow!

Last year an Australian calculator fan sent me improvements on my Romberg-Simpson method. He sent several versions, each slightly better than the one before it. He finally sent me two Excel files that included a good set of test functions. He even included a Timer routine that the VBA code used to measure the time for executing integration routines. I will soon be updating that document and post his Excel test files. You can look at his VBA code for the Romberg method. He performed special test to bring in early termination for the iterations when that was possible. The Excel test files allows you to show or hide the results. Pretty cool stuff.

Namir

[Namir]

I posted the update for the Romberg method. This update includes Excel VBA code and testing from David Graeme.

Namir