Commit a750b6df authored by Lars Stietz's avatar Lars Stietz
Browse files

implementing res and sweep functions

parent ed676b83
......@@ -6,14 +6,31 @@ import time
# Creating quadrature nodes and collocation
num_nodes = 30
num_nodes = 10
t_interval = [0,1]
coll = CollGaussLobatto(num_nodes,t_interval[0],t_interval[1])
# initialize Vector
u_0 = 1
l = -1
u = u_0*np.ones(num_nodes, dtype = np.complex64)
strt = True
nd = True
# if start node is not interval start add value
if coll.nodes[0] != t_interval[0]:
strt = False
u = np.append(u,1)
# if last node is not interval end add value
if coll.nodes[-1] != t_interval[1]:
nd = False
u = np.append(u,1)
# intial data
u_0 = 1
l = -20j
u = np.ones(num_nodes+1, dtype = np.complex64)
l = -1
# RHS of ODE
def f(u):
......@@ -31,6 +48,18 @@ def asmSlowMat(dt):
Q = Q + np.diag(dt[i]*np.ones(sz-i),-i)
return Q
def sweep(u,d_tau,Smat,Qmat,f):
num_nodes = len(u)
u_old = deepcopy(u)
for m in range(1,num_nodes-1):
u[m] = 1/(1-f(1)*d_tau[m])*(u[m-1]-d_tau[m]*f(u_old[m]) + np.dot(Smat[m+1,1:],f(u_old)))
u[num_nodes] = u[0] + np.dot(Qmat[-1,1:],u)
return u
def res(u, Qmat,u_0):
r = u_0 + Qmat[1:,1:]@f(u)-u
return r
solving = 1 #1 solves straight forward, 1 looks at accuracy, 2 compares matrix approach to straight forward
......@@ -38,29 +67,32 @@ if solving == 1:
# SDC straight forward approach
"""
res = np.zeros(num_nodes)
k_max = 20
res_max = np.amax(u[0] +coll.Qmat[1:,1:]@u[1:])
res_max = np.amax(u[0] +coll.Qmat[2 :,:]@u)
res_it = np.empty(0)
j = 0
error = 1
error_it = np.empty(0)
"""
for it in range(250):
#while res_max>1e-15 or error>1e-15:
u_old = deepcopy(u)
for i in range(num_nodes):
u[i+1] = 1/(1-l*coll.delta_m[i])*(u[i]-coll.delta_m[i]*f(u_old[i+1])+np.dot(coll.Smat[i+1][1:],f(u_old[1:])))
for i in range(num_nodes):
res[i] = u[0].real + np.dot(coll.Qmat[i][1:],f(u[1:]).real)-u[i].real
sweep(u,coll.delta_m,coll.Smat,coll.Qmat,f)
r = res(u.real,coll.Qmat,u_0)
"""
for i in range(num_nodes-1):
res[i] = u[0].real + np.dot(coll.Qmat[i][1:-1],f(u[1:-1]).real)-u[i].real
res_max = np.amax(abs(res))
error = np.linalg.norm(u_old.real-u.real)
error_it = np.append(error_it,error)
res_it = np.append(res_it,res_max)
"""
plt.figure(1)
plt.title('SDC solution vs exact solution')
plt.plot(coll.nodes,u[1:].real,'b',label='SDC straight forward')
plt.plot(coll.nodes,u.real,'b',label='SDC straight forward')
plt.plot(coll.nodes,y(coll.nodes).real,'r--',label = 'exact solution')
plt.xlabel('t')
......@@ -69,7 +101,7 @@ if solving == 1:
plt.figure(2)
plt.title('Update ans Residual')
plt.semilogy(error_it,label= 'update difference (2-Norm)')
plt.semilogy(res_it,label='residual')
plt.semilogy(r,label='residual')
plt.xlabel('iteration step')
plt.legend()
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment