Dear All,
When I call the function "G()" which returns two values: p_T and q_T twice but using the same style, that's, P_T, neg = G(); U, Ign = G() and print sums of P_T and U, I get two different answers! Another thing is that both answers are incorrect! Similar problem happens with the function "F()".
Any suggestion of what is wrong and how to fix it?
""" """
from numpy import *
from math import *
from scitools.std import *
N = 1001
M = 200
dx = 2*pi/M
cfl = 0.75
a = 1.0
dt = 0.0005
Eps = 10e-8
U = zeros(M); w = zeros(M); z = zeros(M)
u = zeros(M)
#Solve for the target solution
def Initialise():
return 0.2 + sin(x)
x = linspace(0.0,2*pi, M+1)
U = Initialise()
def Evolve_in_One_Timestep (U):
#create an array to store the new values of the solution
U_new = zeros((M,), float32)
for k in range(1, M-1):
U_new[k] = U[k] - 0.25*(U[k+1]**2 - U[k-1]**2)*(dt/dx) + 0.5*sqrt(a)*(U[k+1] - 2*U[k] + U[k-1])*(dt/dx)
U_new[0] = U_new[1]
U_new[M-1] = U_new[M-2]
return U_new
target = Initialise()
for n in range(0,N):
target = Evolve_in_One_Timestep(target)
#Forward solution
def u_initial_function (x):
value = sin(2*pi*x/M)
return value
u = fromfunction(u_initial_function,(M,))
def f(u):
return u
v = 0.5*(f(u))**2
def Initialize(v):
return v
def Evolve_vs(u,v):
vs = zeros((M,), float32)
for k in range (1, M-1):
vs[k] = Eps/(Eps -dt)*(v[k] - (dt/Eps)*0.5*(u[k])**2)
vs[0] = vs[1]
vs[M-1] =vs[M-2]
return vs
def Evolve_u1(u,v):
u1 = zeros((M,), float32)
vs = Evolve_vs(u,v)
for k in range (1, M-1):
u1[k] = u[k]-0.5*(dt/dx)*(vs[k+1] - vs[k-1]) + 0.5*(dt/dx)*sqrt(a)*(u[k+1] - 2*u[k] + u[k-1])
u1[0] = u1[1]
u1[M-1] = u1[M-2]
return u1
def Evolve_v1(u,v):
v1 = zeros((M,), float32)
vs = Evolve_vs(u,v)
for k in range (1, M-1):
v1[k] = vs[k] - 0.5*(dt/dx)*a*(u[k+1] - u[k-1]) + 0.5*(dt/dx)*(a/sqrt(a))*(vs[k+1] - 2*vs[k] + vs[k-1])
v1[0] = v1[1]
v1[M-1] = v1[M-2]
return v1
def F():
global v ; global u
for n in xrange(N):
u = f(Evolve_u1(u,v))
v = Initialize(Evolve_v1(u,v))
return u, v
#Solve backward problem
p_T , q_T = F()
def V_T(q_T):
return q_T
def U_T(p_T):
return p_T
def EvolveV_T(q_T):
q_ns = zeros((M,), float32)
for k in xrange (M-2, 0,-1):
q_ns[k] = (Eps/(Eps+dt))*q_T[k] + (Eps/(Eps+dt))*(dt/(2.*dx))*(p_T[k+1] - p_T[k-1]) + (Eps/(Eps+dt))*(dt*a/(2.*dx))*(q_T[k+1] - 2.*q_T[k] + q_T[k-1])
q_ns[0] = q_ns[1]
q_ns[M-1] = q_ns[M-2]
return q_ns
def EvolveU_T(p_T):
q_s = EvolveV_T(q_T)
p_ns = zeros((M,), float32)
for k in range (M-2, 0,-1):
p_ns[k] = p_T[k] + (dt/Eps)*q_s[k]*(p_T[k]) + 0.5*(a**2)*(dt/dx)*(q_T[k+1] - q_T[k-1]) + 0.5*a*(dt/dx)*(p_T[k+1] - 2.*p_T[k] + p_T[k-1])
p_ns[0] = p_ns[1]
p_ns[M-1] = p_ns[M-2]
return p_ns
def G():
global p_T; global q_T
for n in range(N):
p_T = U_T(EvolveU_T(p_T))
q_T = V_T(EvolveV_T(q_T))
return p_T, q_T
####################
P_T, neg = G()
print sum(P_T)
U, Ign = G()
print sum(U)
def gradient():
I = sum(dx*(abs(U - target))*P_T)
return I