How to do finite difference integration in python?

Question

I want to solve a differential equation of a given function S. First, I integrate a function using 4th-order runge kutta:

def rk4_step(r, u, h, f, *args):
    k1 = h * f(r, u, *args)
    k2 = h * f(r + h/2, u + k1/2, *args)
    k3 = h * f(r + h/2, u + k2/2, *args)
    k4 = h * f(r + h, u + k3, *args)
    u_new = u + (k1 + 2.0 * k2 + 2.0 * k3 + k4)/6.0
    return u_new

Then I used it to calculate integrate a function and find the S[i] solutions:

S[0] = 0
S[1] = 0
for i in range(1, nr_int):
     S[i+1] = rk4_step(r[i], S[i], dr, rhs_H, index, l, 0, 0, 0, H[i], lamb[i], m[i], nu[ i], p[i], rho[i], cs2[i])
for i in range(nr_int, nr - 1):
     S[i+1] = rk4s_step(r[i], S[i], dr, rhs_H, index, l, 0, 0, 0, 0, lamb[i], m[i], nu[i], 0, 0, 1)

The result is satisfactory and the solutions are coherent.

As a second option, I tried to calculate with finite differences. Being

S[i] = ( S[i+1] - S[i-1] ) / ( 2 * dr )

I used

S[0] = 0
S[1] = 0
for i in range(1, nr_int):
     S[i+1] = S[i-1] + 2 * dr * rhs_H(r[i], S[i], index, l, 0, 0, 0, H[i], lamb[i], m[i], nu[i], p[i], rho[i], cs2[i])
for i in range(nr_int, nr - 1):
     S[i+1] = S[i-1] + 2 * dr * rhs_H(r[i], S[i], index, l, 0, 0, 0, 0, lamb[i], m[i ], nu[i], 0, 0, 1)

But the result printed is very wrong. About 300 orders of magnitude more compared to the calculated method with rouge kutta. What could be wrong with this second method with finite deferences?

300 orders of magnitude is a factor of 10^300. If you actually observe that, then you have divergence. In almost all cases one should stop and re-evaluate the parameters of the integration if the state components change by more than 4 orders of magnitude at the most. (And yes, there is nothing to confirm that Martin Wilhelm K. had a preference for rouge.) — Lutz Lehmann, Oct 23 '22 at 13:35

score 1 · Answer 1 · answered Oct 22 '22 at 15:46

The second method is also known as the central Euler method. It looks nice, symmetric and all, but in the end it is only weakly stable (the stability region is a segment on the imaginary axis), in many cases it has an oscillating error component that grows very fast.

Also your initialization is most likely inappropriate, priming that oscillating error component with a first-order error that gets magnified in the integration steps. It would probably a little more stable if you start the RK4 integration from the index i=0 and use the first value of that in the central Euler integration.

See the section on numerical experiments in https://math.stackexchange.com/questions/3067795/deriving-the-central-euler-method-and-intuition, the relative sizes of the errors for different error orders in the initialization. Note that the error in forward Euler still has order 2.

How to do finite difference integration in python?

1 Answers1