Untitled
unknown
python
3 years ago
9.4 kB
3
Indexable
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator
from numpy.core.function_base import linspace
import scipy.sparse as ssp
from scipy.sparse.linalg import spsolve
from scipy.linalg import circulant
def matrix(N):
dx = 1/(N+1)
matrix = (1/dx**2)*(np.diag(np.ones(N-1), -1) + np.diag(-2*np.ones(N)) + np.diag(np.ones(N-1), 1))
return matrix
def eulerstep(Tdx, uold, dt):
unew = uold + dt*np.dot(Tdx,uold)
return unew
def trapztep(Tdx, uold, dt):
unew = np.linalg.solve(np.eye(len(Tdx)) - np.multiply(Tdx, dt/2), uold + dt/2*np.dot(Tdx,uold))
return unew
def LW_Matrix(N, amu):
# row = np.zeros(shape=(1,N))
# row[0][0] = 1- amu**2
# row[0][1] = (amu/2)*(amu+1)
# row[0][-1] = (amu/2)*(amu-1)
# return ssp.csc_matrix(circulant(row)).todense()
Tdx = (np.diag((amu/2)*(1+amu)*np.ones(N-1), -1) + np.diag((1-np.power(amu,2))*np.ones(N)) + np.diag((amu/2)*(amu-1)*np.ones(N-1), 1))
Tdx[0][-1] = (amu/2)*(1+amu)
Tdx[-1][0] = (amu/2)*(amu-1)
return Tdx
def BTdx(d, dx, sparse ): # without coefficients (dt/(2*dx)) as in own derivation.
N = int(1/dx)
TdxB = (d/dx)*( np.diag(np.ones(N-1),-1) + -2*np.diag(np.ones(N)) + np.diag(np.ones(N-1),1)
)
TdxB[0,-1] = TdxB[1,0]
TdxB[-1,0] = TdxB[0,1]
return TdxB if sparse == False else ssp.csr_matrix(TdxB)
def BurgerMatrix(dx, dt, TdxB, sparse):
N = int(1/dx)
sub = dt - 1
main = (2*dx/dt) - (dt/dx)
sup = dt + 1
Tdx_LW = (dt/(2*dx))*( sub* np.diag(np.ones(N-1),-1)
+ main*np.diag(np.ones(N))
+ sup*np.diag(np.ones(N-1),1)
)
Tdx_LW += (dt/(2*dx))*TdxB
Tdx_LW[0,-1] = Tdx_LW[1,0] # VÄND PÅ MINUSTECKEN OM FEL
Tdx_LW[-1,0] = Tdx_LW[0,1]
return Tdx_LW if sparse == False else ssp.csr_matrix(Tdx_LW)
def BurgerStep(uold, TdxB, burger, d, dx, dt, N):
TdxB[0,-1] = TdxB[1,0]
TdxB[-1,0] = TdxB[0,1]
if(ssp.issparse(TdxB) and ssp.issparse(burger)):
unew = spsolve( (ssp.eye(N) -(d*dt/(2*dx))*TdxB ) , np.dot(burger, ssp.csr_matrix(uold) ) ).toarray()
else:
unew = np.linalg.solve( np.eye(N) - np.multiply(TdxB, (d*dt)/(2*dx)) , np.dot(burger, uold) )
return unew
def LaxWen(u, M):
unew = np.dot(M, u)
return unew
def advectionsSolver(N , M, g, t0, tf):
dx = 1/N
dt = (tf-t0)/M
# a = dx/dt
a = 0.9*dx/dt
amu = a*dt/dx
uold = [g(i/N, i/M) for i in range(1, N+1)]
m = []
m.append(np.array([0, *uold, 0]))
xx = np.linspace(0,1,N+2)
tt = np.linspace(t0,tf,M+1)
T,X = np.meshgrid(tt, xx)
lax = LW_Matrix(N, amu = amu)
for _ in range(1, len(tt)):
uold = LaxWen(uold, lax)
m.append(np.array([0, *uold, 0]))
m = np.array(m)
return m, T, X
def L2(m, tend, taskParam):
errV = []
N = len(m)
t = np.linspace(0,tend,N)
for i in range(0,N):
errV.append(np.sqrt(sum(k*k/(N) for k in m[i])))
plt.plot(t, errV)
plt.xlabel('$dt$')
plt.ylabel('$err$')
plt.title('Error for Custom g(x), with $a \mu = $ ' + taskParam)
plt.show()
def step(choice, Tdx, uold, dt):
if choice == 0:
return eulerstep(Tdx, uold, dt)
if choice == 1:
return trapztep(Tdx, uold, dt)
def plot(m, T, X, task, **param):
fig = plt.figure()
ax = plt.axes(projection = '3d')
wireframe = ax.plot_surface(T,X, np.transpose(m), cmap = cm.plasma) #nipy_spectral, coolwarm , inferno, ocean, jet, gnu, viridis
if task == 1:
title = 'Solution to Diffusion Equation Custom $u(0,x) = g(x) = $'
for i in param:
title = title + param[i] + ' , '
ax.set_title(title)
ax.set_xlabel('T')
ax.set_ylabel('X')
ax.set_zlabel('u(t,x)')
if task == 2:
title = 'Solution to Advection Equation Custom $u(0,x) = g(x,t)$,' +'$a \mu = $'
for i in param:
title = title + param[i] + ' , '
ax.set_title(title)
ax.set_xlabel('T')
ax.set_ylabel('X')
ax.set_zlabel('u(t,x)')
if task == 3:
title = 'Convection Diffusion Custom $u(0,x) = g(x,t)$, '
for i in param:
title = title + param[i] + ' , '
ax.set_title( title )
ax.set_xlabel('T')
ax.set_ylabel('X')
ax.set_zlabel('u(t,x)')
if task == 4:
title = 'Viscous Burger´s Equation, '
for i in param:
title = title + param[i] + ' , '
ax.set_title( title )
ax.set_xlabel('T')
ax.set_ylabel('X')
ax.set_zlabel('u(t,x)')
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter('{x:.02f}')
fig.colorbar(wireframe, shrink = 0.5, aspect = 5)
plt.show()
def task1(N , M, g, t0, tf, choice):
# Next line should compute a 2-dimensional vector u(m∆t,n∆x)
# u(t,x) has u(t,0) = u(t,1) = 0
# u(0,x) = g(x)
dt = (tf-t0)/M
xx = np.linspace(0,1,N+2)
tt = np.linspace(t0,tf,M+1)
T,X = np.meshgrid(tt, xx)
uold = [g(i/(N), i/(M)) for i in range(1, N+1)]
m = []
m.append(np.array([0, *uold, 0]))
for _ in range(1, np.size(tt)):
uold = step(choice, Tdx = matrix(N), uold = uold, dt = dt )
m.append(np.array([0, *uold, 0]))
m = np.array(m)
return m, T, X
def convdif(a,d,dt,dx): # u' = (d*Tdx - a*Sdx)u
N = int(1/dx)
M = int(1/dt)
g = lambda x,t: x**3 -(4/3)*x**2 +(1/3)*x
uold = [g(i/N, i/M) for i in range(1, N+1)]
xx = np.linspace(0,1,N)
tt = np.linspace(0,1,M+1)
T,X = np.meshgrid(tt, xx)
Tdx = (np.diag(1*np.ones(N-1), -1) + np.diag(-2*np.ones(N)) + np.diag(1*np.ones(N-1), 1))/dx**2
Sdx = (np.diag(-1*np.ones(N-1), -1) + np.diag(0*np.ones(N)) + np.diag(np.ones(N-1), 1))/(2*dx)
Master = (d*Tdx - a*Sdx)
Master[0][-1] = d/(dx**2)+a/(2*dx)
Master[-1][0] = d/(dx**2)-a/(2*dx)
m = []
m.append(np.array([*uold]))
for _ in range(0, np.size(tt)-1):
uold = step(1, Master, uold, dt)
m.append(np.array([*uold]))
m = np.array(m)
print(np.shape(np.transpose(m)))
return m, T, X
def viscBurgerSolver(g, d):
N = 300
delx = 1/N
M = 5*N
delt = 1/M
uold = [g(i/N, i/M) for i in range(1, N+1)]
BTedX = BTdx(d,delx, sparse = False )
Burger = BurgerMatrix(delx, delt, BTedX, sparse = False)
m = []
m.append(np.array([*uold]))
for _ in range(0, M):
uold = BurgerStep(uold,BTedX, Burger, d, delx, delt, N)
m.append(np.array([*uold]))
m = np.array(m)
xx = np.linspace(0,1,N)
tt = np.linspace(0,1,M+1)
T,X = np.meshgrid(tt, xx)
print(np.shape(m))
return m, T , X
if __name__ == '__main__':
c = 1e-6
N = 200
tstart = 0
tend= 1
t = tend - tstart
M_CFL = 2*t*(N+1)**2 # CFL-condition dt/dx**2 <= 1/2 \iff t/M * (N+1)^2
chc = 1
if chc == 0:
Mnbr = int(np.ceil(M_CFL))
if chc == 1:
Mnbr = 10*N
# m, T, X = task1(N = 100, M = Mnbr, g = lambda x,t: (x-0.5)**2, t0= tstart, tf=tend, choice = chc )
# m, T, X = task1(N = N, M = Mnbr-10, g = lambda x,t: c*(x)**3 + c/((t)**3), t0= tstart, tf=tend, choice = chc) #BREAK IT!
# m, T, X = task1(N = N, M = Mnbr, g = lambda x,t: np.exp(x /t), t0= tstart, tf = tend, choice = chc )
# m, T, X = task1(N = N, M = Mnbr-16, g = lambda x,t: np.exp(x /t), t0= tstart, tf = tend, choice = chc) #BREAK IT AND GET PANTS!
# plot(m = m, T = T, X = X, task=1)
# advectionsSolver(N = N, M = Mnbr, g = lambda x: np.exp(-100*(x-0.5)**2))
# m, T, X = advectionsSolver(N = 100, M = Mnbr, g = lambda x,t: np.exp(-100*np.power((x-0.5),2)), t0 = 0, tf = 1 )
# plot(m = m, T = T, X = X)
#a = -10
#d = a/(N**2)
# d = 1000000
# print('Pe') # d > a * dx/2 = a ^(1/N^2),
# print(np.abs(d/a))
# m, T, X = convdif(a, d, 1/Mnbr, 1/(N+1) )
# plot(m, T, X, task = 3, a_and_d = '$(a,d) = (1, 1.11e-5)$', g_of_x = '$g(x) = x^3 -x^2 4/3 + x/3 $')
# L2(m,5, taskParam='$0.9$')
##########UPPGIFT 4.1 / 4.1###########################
# Jag tror mina matriser är fel, har tolkat LW(uold) som en matrismultiplikation U*uold
# där U härleds till (dt/2dx)*tridiag([dt-1, (2*dx/dt) - (dt/dx) ,dt+1 ])
# + de hörnens punkter så det blir en circulent matris. Jag har hittat att U:s invers är exakt U^T fast alla nollskilda element byts ut till 1/a_ij. Inte till någon hjälp precis då det är (I - TdxCirc)^-1 som ska lösas i BurgerStep.
# Iförsig är den matrisen cirkulant med
# Circ( 1/(1 + 2/dx**2), -dx**2, 0, ... , 0, -dx**2 )!!!!!
# Vi testar imorgon godnatt!!
# Sparse fungerar ej... Men vi behöver omimplementera matriserna och kanske lösa den bättre
m, T, X = viscBurgerSolver(lambda x,t: np.exp(-100*np.power((x-0.5),2) ), d = 0.01) #Lägre d fungerar men inte 0.1..
plot(m, T, X, 4, dcoeff = '$d = 0.01$', g_of_x = '$g(x) = \exp(-100 (x-0.5)^2)$' )
Editor is loading...