# Project1

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()

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()

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()

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__':