Мой друг и я написали следующий код для решения уравнений с частными производными:
import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sp
import scipy.sparse.linalg as la
import gc
import time
#%%
h = 0.1
const = 1
dt1 = 0.25*0.99*h**2
dt2 = 0.1
dt3 = 0.001
Nt_1 = int(40.0/dt1+1)
Nt_2 = int(40.0/dt2+1)
Nt_3 = int(40.0/dt3+1)
t_1 = np.linspace(0,40.0,Nt_1)
t_2 = np.linspace(0,40.0,Nt_2)
t_3 = np.linspace(0,40.0,Nt_3)
xbegin = 0
xend = 16
ybegin = 0
yend = 8
Sx = xend-xbegin
Sy = yend-xbegin
Nx = int(Sx/h+1)
Ny = int(Sy/h+1)
Np = (Nx-2)*(Ny-2)
I = sp.eye(Np,Np,format='csc')
def meshcreator(xb,xe,yb,ye):
meshy1,meshx1 = np.mgrid[yb+h:ye:h,xb+h:xe:h]
return meshy1,meshx1
def Matrixmaker(Nx3,Ny3):
Dx = 1/h*(sp.diags([-1, 1], [-1, 0], shape=(Nx3-1, Nx3-2),format='csc'))
Dy = 1/h*(sp.diags([-1, 1], [-1, 0], shape=(Ny3-1, Ny3-2),format='csc'))
Dx_t = Dx.transpose()
Dy_t = Dy.transpose()
Lxx = Dx_t.dot(Dx)
Lyy = Dy_t.dot(Dy)
Ix = sp.eye(Nx3-2)
Iy = sp.eye(Ny3-2)
L = sp.kron(Iy,Lxx) + sp.kron(Lyy,Ix)
L = -L
return L
#%%
def k_maker(h1,Nx1,Ny1,q):
meshy2,meshx2 = meshcreator(xbegin,xend,ybegin,yend)
# k1 = sp.lil_matrix((Ny1-2, Nx1-2)).toarray()
k1 = np.zeros((Ny1-2,Nx1-2))
for i in range(len(meshx2[:,0])):
for j in range(len(meshx2[0,:])):
if 1 <= meshx2[i,j] <= 2 and 1 <= meshy2[i,j] <= 2:
k1[i,j] = q
elif 1 <= meshx2[i,j] <= 3 and 3 <= meshy2[i,j] <= 5:
k1[i,j] = q
elif 4 <= meshx2[i,j] <= 7 and 4 <= meshy2[i,j] <= 7:
k1[i,j] = q
elif 9 <= meshx2[i,j] <= 12 and 4 <= meshy2[i,j] <= 6:
k1[i,j] = q
elif 13 <= meshx2[i,j] <= 15 and 1 <= meshy2[i,j] <= 3:
k1[i,j] = q
return k1, meshx2, meshy2
k,meshx,meshy = k_maker(h,Nx,Ny,const)
plt.figure(4)
plt.spy(k,aspect='equal')
k = k.reshape((Np,1))
#%%
def ut0(x,y):
return (np.exp(-2*(x-1.5)**2-2*(y-1.5)**2).reshape(Np))
def jacobian(L2,value):
value2 = k*value
u_diag = sp.diags(diagonals = [np.reshape(value2,(Np))],offsets=[0],shape=(Np,Np),format='csc')
jacob = L2 + sp.diags(diagonals = [k.reshape(Np)],offsets=[0],shape=(Np,Np),format='csc') - 2*u_diag
return jacob
#%%
A = Matrixmaker(Nx,Ny)
ut0 = (ut0(meshx,meshy).reshape(Np,1))
def FE_BE_iterator(L3,u1,u0,timestep,Nt,choice,t,eps1 = 10e-3):
def f1(uk):
return L3*uk
def f2(uk):
return k*uk
def f3(uk):
utmp = uk*uk
return k*utmp
def F0(uk):
return f1(uk)+f2(uk)-f3(uk)
errorlist = []
usollist = [u0]
count = 0
if choice == 'Forward':
timerstart = time.time()
for i in range(Nt):
usol = timestep*(F0(u0))+u0
u0 = usol
if abs(t[i]-1) <= timestep and count == 0:
usollist.append(usol.reshape(Np,1))
count +=1
if abs(t[i]-2) <= timestep and count == 1:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-3) <= timestep and count == 2:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-5) <= timestep and count == 3:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-10) <= timestep and count == 4:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-20) <= timestep and count == 5:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-40) <= timestep and count == 6:
usollist.append(usol.reshape(Np,1))
count+=1
timerend = time.time()
time_elapsed = timerend-timerstart
print("Time elapsed for " + choice + " = ", time_elapsed)
return usollist,errorlist
if choice == 'Newton':
timerstart = time.time()
for i in range(Nt):
residual = 1
while residual > eps1:
jac = jacobian(L3,u1)
S = I-jac*timestep
RHS = u1-u0-timestep*F0(u1)
v = (la.spsolve(S,RHS)).reshape(Np,1)
ukp1 = u1-v
u1 = ukp1
residual = np.linalg.norm(v)
errorlist.append(residual)
usol = timestep*(F0(ukp1)) + u0
u0 = usol
if abs(t[i]-1) <= timestep and count == 0:
usollist.append(usol.reshape(Np,1))
count +=1
if abs(t[i]-2) <= timestep and count == 1:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-3) <= timestep and count == 2:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-5) <= timestep and count == 3:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-10) <= timestep and count == 4:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-20) <= timestep and count == 5:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-40) <= timestep and count == 6:
usollist.append(usol.reshape(Np,1))
count+=1
timerend = time.time()
time_elapsed = timerend-timerstart
print("Time elapsed for " + choice + " = ", time_elapsed)
return usollist,errorlist
if choice == 'Picard':
timerstart = time.time()
for i in range(Nt):
usol = u0
residual = 1
while residual > eps1:
usol = timestep*F0(usol) + u0
residual = np.linalg.norm(usol-timestep*F0(usol)-u0)
u0 = usol
if abs(t[i]-1) <= timestep and count == 0:
usollist.append(usol.reshape(Np,1))
count +=1
if abs(t[i]-2) <= timestep and count == 1:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-3) <= timestep and count == 2:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-5) <= timestep and count == 3:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-10) <= timestep and count == 4:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-20) <= timestep and count == 5:
usollist.append(usol.reshape(Np,1))
count+=1
if abs(t[i]-40) <= timestep and count == 6:
usollist.append(usol.reshape(Np,1))
count+=1
timerend = time.time()
time_elapsed = timerend-timerstart
print("Time elapsed for " + choice + " = ", time_elapsed)
return usollist,errorlist
def plot(M,start,startp1,stepdt,N_t,name,t_time):
Solution1 = FE_BE_iterator(M,start,startp1,stepdt,N_t,name,t_time)
Solution = Solution1[0]
Error = Solution1[1]
for s in range(len(Solution)):
plt.subplot(2,4,s+1)
plt.imshow((Solution[s]).reshape((Nx-2,Ny-2),order ='F'),origin='lower',extent=(ybegin,yend,xbegin,xend))
plt.colorbar()
if name == 'Newton':
plt.figure(0)
plt.plot(Error)
counter = 0
for j in ['Forward','Newton','Picard']:
plt.figure(counter+1)
if j == 'Forward':
plot(A,ut0,ut0,dt1,Nt_1,j,t_1)
elif j == 'Newton':
plot(A,ut0,ut0,dt2,Nt_2,j,t_2)
elif j == 'Picard':
plot(A,ut0,ut0,dt3,Nt_3,j,t_3)
counter+=1
gc.collect()
Это дает такие графики, как этот: 
Принимая во внимание, чтоЯ хотел бы, чтобы они выглядели так:
Я знаю, что вектор решения в функции imshow должен быть транспонирован, но транспонировать егопросто портит мое решение и не поворачивает изображение. Как это исправить?
РЕДАКТИРОВАТЬ: Вот что дает моя попытка поворота изображения с помощью транспонирования: 