-
Notifications
You must be signed in to change notification settings - Fork 0
/
L2ConvergenceAdvection.py
72 lines (59 loc) · 1.75 KB
/
L2ConvergenceAdvection.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
from math import sin,sqrt
import numpy as np
import finitedifferences
from matplotlib import pyplot as plt
import time
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from mpl_toolkits import mplot3d
PI = 3.1415
def uExc(u,T,k,x):
uExc = u.copy()
for i in range(1,u.shape[1]):
t=k*i
uExc[:,i]=np.sin(2*3.1415*(x + t))
return uExc
def RK4(A,f,u,h,k):
u[:,0]=f[:]
for i in range(1,u.shape[1]):
k1 = k/h*A @ u[:,i-1]
k2 = k/h*A @ (u[:,i-1] + 0.5*k1)
k3 = k/h*A @ (u[:,i-1] + 0.5*k2)
k4 = k/h*A @ (u[:,i-1] + k3)
u[:,i] = u[:,i-1] + 1/6*(k1 + 2*k2 + 2*k3 + k4)
return u
def L2Error(m,n,T):
u = np.zeros([m,n])
h = 1/m
k = T/(n-1)
A = finitedifferences.centralDifference(m)
x = np.linspace(0,1,m,endpoint= False)
f = np.sin(2*PI*x)
u = RK4(A,f,u,h,k)
#uE = uExc(u,T,k,x)
error = np.abs(u[:,n-1]-f)
return sqrt(np.dot(error,error)*h)
nTests = 7
gridSize = 4
res = np.zeros([nTests,4])
for i in range(0,nTests):
res[i,0] = L2Error(gridSize,5*2*2*2**nTests,1)
res[i,1] = 2/gridSize
if(i>0):
[res[i,2], intercept] = np.polyfit(np.log(res[i-1:i+1,1]),np.log(res[i-1:i+1,0]),1)
gridSize = gridSize * 2
[slope, intercept] = np.polyfit(np.log(res[:,1]),np.log(res[:,0]),1)
print(res)
plt.figure(num=None, figsize=(12, 6), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(121)
plt.plot(res[:,1],res[:,0],"bs")
plt.xlabel('steplength in space and time')
plt.ylabel('L2 error')
plt.gca().invert_xaxis()
plt.subplot(122)
plt.plot(res[1:,1],res[1:,2],"g^")
plt.xlabel('steplength in space and time')
plt.ylabel('Convergence/slope of loglog plot')
plt.gca().invert_xaxis()
plt.savefig("loglog_slope.png")
plt.show()