0
$\begingroup$

I'm trying to estimate pi by solving the IVP
$y'' + y = 0$
where
$y(0) = 1, y'(0) = 0$
numerically by defining $\frac{\pi}{2}$ as the first value on t such that $y(t) = 0$

I'm trying to solve this problem both with Euler forward method and with Runge-Kutta 4 method.
However as of right now I get the exact same value from both Euler and RK4 which I think is weird since the error should be higher in the Euler method. (?)

Right now I'm stuck since I don't know what's causing this. I think there could be either one of two problems (maybe both):

  1. My RK4 is somehow wrong
  2. The way I try to estimate pi is wrong

I can't find any problem with my RK4 method but I'm not overly confident with this method.. yet.
I think (2) could be the error but if that's the case I don't know how I'm supposed to estimate pi.
I'm aware of the Newton-Raphson method and the bisection method but both of them require a function.. all I have is y values & t values, not a function.
Right now all I'm doing is stopping the simulation as soon as $y(t)$ drops below zero at which I assume t is close enough to $\pi$, which is the first root.

If anyone has any ideas regarding what I'm doing wrong or how I could fix this I'd be very happy.

I'll paste my code (in python) below:

import numpy as np

'''
D.E:
y' = u
u' = -y
'''
def f(u):
    return np.array([u[1],-u[0]])

# Euler:
def Euler(y0, t, dt):
    t0 = t[0]
    for _ in range(1, len(t)):
        y0 = y0+dt*f(y0) 
        if y0[0]<0: # I assume that we're close to pi if y drops below zero and if that's the case we want to exit the loop
            return y0,t0
        t0 = t0 + dt
    return y0, t0

# RK 4
def RK4step(u,dt):
    k1 = f(u)
    k2 = f(u+0.5*k1*dt)
    k3 = f(u+0.5*k2*dt)
    k4 = f(u+k3*dt)
    return u + dt*(k1+2*k2+2*k3+k4)/6

def RK4(y0, t, dt):
    t0 = t[0]
    for _ in range(1, len(t)):
        y0 = RK4step(y0,dt)
        if y0[0] < 0: # I assume that we're close to pi if y drops below zero and if that's the case we want to exit the loop
            return y0,t0
        t0 = t0 + dt
    return y0, t0

if __name__ == "__main__":
    real_pi = np.pi

    # Initial values
    T1 = np.pi/2
    T2 = 4
    u0 = np.array([1, 0]) # Initial values
    n = 5

    # Want to loop for different lengths of dt, ideally more than 5 but the program gets very slow after 5.
    for i in range(n):
        dt = 10**-(i+1)
        n = int((T2-T1)/dt)
        t = np.linspace(T1,T2,n)

        y_euler,t_euler = Euler(u0, t, dt)
        y_RK4,t_RK4 = RK4(u0, t, dt)

        print()
        print(f"{i+1}:")
        print(f"------Euler:------")
        print(f"Estimated pi: {t_euler}")
        print(f"Error Euler = {np.abs(t_euler-real_pi)}")

        print(f"------Runge-Kutta 4:------")
        print(f"Estimated pi: {t_RK4}")
        print(f"Error RK4 = {np.abs(t_RK4-real_pi)}")
```
$\endgroup$

1 Answer 1

1
$\begingroup$

Set T1=0, you want to approximate $\pi$ without knowing it. Just multiply the result with 2.

This will still give the same result in both methods. The reason is that the angle in the Euler method has a quadratic error, each step has an angle of $\arctan(h)=h+\frac{h^3}3+...$, so the accumulation of $n=t/h$ steps gives $n\arctan(h)=t+\frac{th^2}3+...$ The RK4 step is of course even more exact. Thus both methods pass the zero line after the same number of steps, except in some exceptional cases that were not encountered.

To get a better comparison add an estimate of where in the step the root is. Use a Newton step to refine the intersection point, using that the state vector is $[y,y']$. Using the previous point one could also use the secant root or a reverse cubic interpolation.

Perform the state and time updates at the same time to remain consistent. The loop would then look like

def odeint(stepper, y0, t, dt):
    t0 = t[0]
    for _ in range(1, len(t)):
        y0, t0 = stepper(y0,dt), t0 + dt
        if y0[0] < 0: # I assume that we're close to pi if y drops below zero and if that's the case we want to exit the loop
            break
    return y0, t0 - y0[0]/y0[1]

...

        y_euler,t_euler = odeint(Eulerstep, u0, t, dt)
        y_RK4,t_RK4 = odeint(RK4step, u0, t, dt)

With that RK4 reaches full accuracy of an relative error of $10^{-12}$ in the third step

1:
------Euler:------
Estimated pi: 3.1521866680981585
Error Euler = 0.01059401450836539
------Runge-Kutta 4:------
Estimated pi: 3.141578703016505
Error RK4 = 1.395057328812399e-05

2:
------Euler:------
Estimated pi: 3.1416974697101296
Error Euler = 0.0001048161203365261
------Runge-Kutta 4:------
Estimated pi: 3.14159213408815
Error RK4 = 5.195016430370458e-07

3:
------Euler:------
Estimated pi: 3.141593700916784
Error Euler = 1.0473269909105909e-06
------Runge-Kutta 4:------
Estimated pi: 3.1415926535840635
Error RK4 = 5.729638985485508e-12
$\endgroup$
1
  • $\begingroup$ Thanks a lot for the help! I appreciate it. $\endgroup$
    – Minea
    Oct 10, 2022 at 20:15

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .