Newton's method, Trapezoidal method, Euler's method, Runge-Kutta 2, Runge-Kutta 4

numericalapproximation.py

import numpy as np
from matplotlib import pyplot as plt

def f(l, a, b, T):
    return a*(l**2) - a/l + b*l - b/(l**2) - T

def df(l, a, b):
    return 2*a*l + a/(l**2) + b + 2*b/(l**3)

def q1_newton(x0, f, df, a, b, T, n=1):
    print(f"x0 = {x0}.")
    x = x0
    for i in range(n):
        x = x-f(x, a, b, T)/df(x, a, b)
        print(f"x{i+1} = {x}")
    return x

def q1_trapezoid(l, a, b, n):
    d = (b-a)/n
    T = []
    for i in range(n):
        T.append(d/2*(l[i+1]+l[i]))
    return sum(T)

def T_compute(x, a, b, n):
    T = []
    for i in range(n):
        T.append((x[i]**2) * (a*(1+(x[i]**2))+b) / (1+(x[i]**2))**2)
    return T

def u_compute(x,f,df,a,b,l,n,u0):
    x0 = 1
    T = T_compute(x,a,b,n)
    l = [x0]
    for i in range(n):
        l.append(q1_newton(x0,f,df,a,b,T[i]))
    u = [u0]
    for i in range(n):
        u.append(u0-x[i+1]+q1_trapezoid(l,0,x[i+1],i+1))
    return u


xlist = np.linspace(0,1,21)
ulist = u_compute(xlist,f,df,a=20,b=10,l=1,n=20,u0=0)

plt.scatter(xlist,ulist)
plt.xlabel("x")
plt.ylabel("u")
plt.title("q1 plot")
plt.show()

print(xlist)
print(ulist)

print("-----------------------------")


def euler(t0,tf,y0,n):
    h = (tf-t0)/(n-1)
    t=np.linspace(t0,tf,n)
    y = np.zeros([n])
    y[0]=y0
    print(f"t[0],y[0] = {t[0]}, {y[0]}")
    for i in range(1,n):
        y[i] = h*(np.exp(-y[i-1]) + t[i-1]**5) + y[i-1]
        print(f"t[{i}],y[{i}] = {t[i]}, {y[i]}")

    plt.plot(t,y,'o')
    plt.xlabel("x")
    plt.ylabel("y_approximated (Euler)")
    plt.title("Q2: Euler's approximations")
    plt.show()


def rk2(t0,tf,y0,n):
    h = (tf-t0)/(n-1)
    t=np.linspace(t0,tf,n)
    y = np.zeros([n])
    y[0]=y0
    print(f"t[0],y[0] = {t[0]}, {y[0]}")
    for i in range(1,n):
        k1 = h * (np.exp(-y[i-1]) + t[i-1]**5)
        k2 = h * (np.exp(-y[i-1]-k1) + (t[i])**5)
        y[i] = y[i-1] + ( k1 + k2 ) / 2.0
        print(f"t[{i}],y[{i}] = {t[i]}, {y[i]}")

    plt.plot(t,y,'o')
    plt.xlabel("x")
    plt.ylabel("y_approximated (Runge-Kutta 2)")
    plt.title("Q2: RK2 approximations")
    plt.show()

def rk4(t0,tf,y0,n):
    h = (tf-t0)/(n-1)
    t=np.linspace(t0,tf,n)
    y = np.zeros([n])
    y[0]=y0
    print(f"t[0],y[0] = {t[0]}, {y[0]}")
    for i in range(1,n):
        k1 = h * (np.exp(-y[i-1]) + (t[i-1])**5)
        k2 = h * (np.exp(-(y[i-1] + 0.5 * k1)) + (t[i-1] + 0.5 * h)**5)
        k3 = h * (np.exp(-(y[i-1] + 0.5 * k2)) + (t[i-1] + 0.5 * h)**5)
        k4 = h * (np.exp(-(y[i-1] + k3)) + (t[i])**5)
        y[i] = y[i-1] + ( k1 + 2.0 * ( k2 + k3 ) + k4 ) / 6.0
        print(f"t[{i}],y[{i}] = {t[i]}, {y[i]}")

    plt.plot(t,y,'o')
    plt.xlabel("x")
    plt.ylabel("y_approximated (Runge-Kutta 4)")
    plt.title("Q2: RK4 approximations")
    plt.show()


t0=0
y0=1
tf=1
n=100

rk4(t0,tf,y0,n)
rk2(t0,tf,y0,n)
euler(t0,tf,y0,n)