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
```