Project1

mail@pastecode.io avatar
unknown
python
2 years ago
6.4 kB
1
Indexable
Never
from functools import total_ordering
import numpy as np
import matplotlib.pyplot as plt
from numpy.lib.twodim_base import vander
from scipy.linalg import expm
from scipy import integrate

a = 3
b = 9
c = 15
d = 15

lamb = 1
def testeq(_,u): 
    l = np.array([[lamb]])
    return l*u

def testeqnbr(_,u): 
    return lamb*u

def solution(y0, tg):
    l = np.array([[lamb]])
    return np.array([np.dot(np.exp(l*t), y0) for t in tg])

def RK34Step(f, told, uold, h):
    unew = RK4Step(f, told, uold, h)
    hY1 = h*f(told,uold)
    hY2 = h*f(told+h/2, uold+hY1/2)
    hY3 = h*f(told+h/2, uold+hY2/2)
    hY4 = h*f(told+h, uold+hY3)
    hZ3 = h* f(told + h, uold - hY1 + 2*hY2)
    err = (2*hY2 + hZ3 -2*hY3 - hY4)/6
    return [unew, err]
    

def RK4Step(f, told, uold, h): 
    hY1 = h*f(told,uold)
    hY2 = h*f(told+h/2, uold+hY1/2)
    hY3 = h*f(told+h/2, uold+hY2/2)
    hY4 = h*f(told+h, uold+hY3)
    return uold + (hY1 + 2*hY2 + 2*hY3 + hY4)/6

def newstep(tol, err, errold, hold, k):
    return ( (tol/abs(err))**(2/(3*k)) * (tol/abs(errold))**(-1/(3*k)) )*hold

def adaptiveRK34(f, t0, tf, y0, tol):
    t = [t0]
    y = [y0]
    
    tn = t0
    yn = y0
    rn = tol
    hm = abs(tf-t0)*np.power(tol, 1/4)/(100*(1 + np.linalg.norm(f(None, y0))))
   
    while(hm < tf - tn):
        yn, err = RK34Step(f, tn, yn, hm)
        tn= tn+hm

        t.append(tn) 
        y.append(yn)
        
        rold = rn
        rn = np.linalg.norm(err)
        hm = newstep(tol, rn, rold, hm, 4)
       

    t.append(tf) 
    y.append(RK34Step(f, tn, yn, tf - tn)[0] )
    return [np.array(t), np.array(y)]

def RK4Solve(f, t0, tf, y0, N):
    tg = np.linspace(t0,tf,N+1)
    h =(tf-t0)/N
    approx=np.array([y0])
    yn = y0
    tn = t0
    for _ in range(N):
        yn = RK4Step(f, tn, yn, h)
        tn = tn+h
        approx = np.append(approx, yn)
   
    err = solution(y0, tg) - approx
    return [tg, approx, err]

def lotkaVolterraSolver():
    u_f0 = np.array([1,1])
    u2_f0 = np.array([1.7, 0.1])
    tol = 1e-7
    tmax = 10
    lotka = lambda t, u: np.array([a*u[0] - b*u[0]*u[1], c*u[0]*u[1] - d*u[1] ])
    [t, y] = adaptiveRK34(lotka, 0, tmax, u_f0, tol)
    [t2, y2] = adaptiveRK34(lotka, 0 , tmax, u2_f0, tol)

    title = "Solution to initial condition over time t with u_f0 = [ "
    add = ",".join(str(x) for x in u_f0)
    add2 = ",".join(str(x) for x in u2_f0)

    figure, axis = plt.subplots(2,2)
    axis[0, 0].plot(t, y)
    axis[0, 0].set_title(title + add + " ]")

    axis[0, 1].plot(y[:,0], y[:,1], 'm')
    axis[0, 1].set_title("Phase portrait when u_f0 = [ " + add + " ]")

    axis[1, 0].plot(t2, y2)
    axis[1, 0].set_title(title + add2 + " ]")

    axis[1, 1].plot(y2[:,0], y2[:,1], 'm')
    axis[1, 1].set_title("Phase portrait when u_f0 = [ " + add2 + " ]")
    plt.show()

    plt.figure(2)
    plt.title("H:s change over time")
    plt.semilogy(t, [abs(H(u ,a ,b ,c ,d)/H(u_f0 , a, b ,c ,d) - 1) for u in y], 'r')
    plt.show()

def H(u, *parameters):
    x, y = u
    a, b, c, d = parameters
    sol = c*x + b*y - d*np.log(x) - a*np.log(y)
    return sol
    
def errVSh(f, t0, tf, y0):
    nbr = 11
    N = 2**np.linspace(1, nbr, num=nbr, dtype=int)
    h = np.divide(tf-t0, N)
    E = []
    for n in N:
        _, _, err = RK4Solve(f,t0,tf,y0,n)
        E.append(np.linalg.norm(err[-1][-1])) #ska vara två [][] för matris!!!!!!!!
    plt.figure(1)
    plt.loglog(h,E)
    plt.loglog(h, h**4)
    plt.title("Error plot")
    plt.show()

def task1Solver():
     y0 = np.array([1])
     t0 = 0
     tf = 1
     [tg, app, err] = RK4Solve(testeq,t0,tf,y0, 100)
     print(err)
     errVSh(testeq, 0 ,1 ,y0)
     t = np.linspace(0, 1 ,1000)
     plt.figure(2)
     plt.plot(tg, app, 'k')    
     plt.plot(t, solution(y0,t), 'r--' )
     plt.show()

     t,y = adaptiveRK34(testeqnbr, 0, 1, 1, 10**(-6))  
     tbig = np.linspace(0,1,1000)  


     plt.plot(t,y)
     plt.plot(tbig, solution(np.array([1]), tbig))
     plt.show()

def VanDerPolSolver():
    mu = 100
    tol = 1e-7
    t0 = 0 
    tf = 300
    u_f0 = np.array([2, 0])
    u2_f0 = np.array([3 , 1]) 
    vanDerPol = lambda t, u: np.array([u[1] , mu*(1-u[0]**2)*u[1] -u[0] ])
    [t, y] = adaptiveRK34(vanDerPol, t0, tf, u_f0, tol)
    [t2, y2] = adaptiveRK34(vanDerPol, t0, 3*tf, u2_f0, tol)

    title = "Solution to initial condition over time t with u_f0 = [ "
    add = ",".join(str(x) for x in u_f0)
    add2 = ",".join(str(x) for x in u2_f0)

    figure, axis = plt.subplots(2,2)
    axis[0, 0].plot(t, y[:,1])
    axis[0, 0].set_title(title + add + " ]")

    axis[0, 1].plot(y[:,0], y[:,1], 'm')
    axis[0, 1].set_title("Phase portrait when u_f0 = [ " + add + " ]")

    axis[1, 0].plot(t2, y2[:,1], 'r')
    axis[1, 0].set_title(title + add2 + " ]")

    axis[1, 1].plot(y2[:,0], y2[:,1], 'm')
    axis[1, 1].set_title("Phase portrait when u_f0 = [ " + add2 + " ]")
    plt.show()
    #Task 3.2
    
def VanDerPol32():
    muus = np.array([10, 15, 22, 33, 47, 68, 100, 150, 220, 330, 470]) 
    tol = 1e-7
    t0 = 0 
    u_f0 = np.array([2, 0])
   
    Ns = []
    for i in range(0,len(muus),1):
        vanDerPol = lambda t, u: np.array([u[1] , muus[i]*(1-u[0]**2)*u[1] -u[0] ])
        t, _,  = adaptiveRK34(vanDerPol, t0, 0.7*muus[i], u_f0, tol)
        N = len(t)-1
        Ns.append(N)
        
    

    plt.figure(1)
    plt.plot(muus,Ns, 'm')
    plt.title("Stiffness: Number of Steps Against Müh")
    plt.xlabel('Müh')
    plt.ylabel('N')
    plt.show()

    #Task 3.3 only
    
    muus2 = np.array([10, 15, 22, 33, 47, 68, 100, 150, 220, 330, 470, 680, 1000, 4000, 10000])
    Ns2 = []
    for i in  range(0,len(muus2),1):
        vanDerPol = lambda t, u: np.array([u[1] , muus2[i]*(1-u[0]**2)*u[1] -u[0] ])
        tgrid = integrate.solve_ivp(vanDerPol, [0, 0.7*muus2[i]],  u_f0, method = 'BDF').t
        N = len(tgrid)-1
        Ns2.append(N)

    plt.figure(2)
    plt.semilogx(muus2, Ns2)
    plt.title("Only implicit BDF Method")
    plt.xlabel('Müh')
    plt.ylabel('N')
    plt.show()


if __name__ == '__main__':
    task1Solver()
    # lotkaVolterraSolver()
    # VanDerPolSolver()
    # VanDerPol32()