Untitled
unknown
python
3 years ago
6.8 kB
5
Indexable
from functools import total_ordering
import numpy as np
import matplotlib.pyplot as plt
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): #f är lambdauttryck motsvarande, tar in 2 parametrar
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))))
nbrSteps = 0
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)
nbrSteps += 1
t.append(tf)
y.append(RK34Step(f, tn, yn, tf - tn)[0] )
nbrSteps += 1
return [np.array(t), np.array(y), nbrSteps]
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 = 100
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, 680, 1000]) #1000 läggs till så blir det 13 fönster
tol = 1e-7
t0 = 0
tf = 50
u_f0 = np.array([2, 0])
Ns = []
Ns2 = []
for i in range (0, 13,1):
vanDerPol = lambda t, u: np.array([u[1] , muus[i]*(1-u[0]**2)*u[1] -u[0] ])
_, _, nbr = adaptiveRK34(vanDerPol, t0, tf, u_f0, tol)
Ns.append(nbr)
sol = integrate.solve_ivp(vanDerPol, [t0,tf], u_f0, method='BDF' )
implicit = len(sol.t) - 1
Ns2.append(implicit)
print(Ns)
print(Ns2)
plt.figure(1)
plt.plot(muus,Ns, 'm')
plt.plot(muus, Ns2, 'b')
plt.title("Stiffness: Number of Steps Against Müh")
plt.show()
#Task 3.3 only
muus2 = np.array([10, 15, 22, 33, 47, 68, 100, 150, 220, 330, 470, 680, 1000, 4000, 10000])
print(muus2)
vanDerPol4000 = lambda t, u: np.array([u[1] , 4000*(1-u[0]**2)*u[1] -u[0] ])
sol4000 = integrate.solve_ivp(vanDerPol4000, [t0,tf], u_f0, method='BDF' )
implicit = len(sol4000.t) - 1
Ns2.append(implicit)
vanDerPol10000 = lambda t, u: np.array([u[1] , 10000*(1-u[0]**2)*u[1] -u[0] ])
sol10000 = integrate.solve_ivp(vanDerPol10000, [t0,tf], u_f0, method='BDF' )
implicit = len(sol10000.t) - 1
Ns2.append(implicit)
plt.figure(2)
plt.semilogy(muus2, Ns2)
plt.title("Only implicit BDF Method")
plt.show()
if __name__ == '__main__':
# task1Solver()
# lotkaVolterraSolver()
# VanDerPolSolver()
VanDerPol32()
Editor is loading...