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

Clean up your code, if you want response, by using sensible variable names, not using meaningless functions like

def Initialize(v):
	return v

which is equivalent of simple v and does not create Initialize object like name suggest.

Also please put first imports, global constants, then function definitions and at the end the main routine either as one more function called from main level or simply at main level, but do not intersperse function definitions with calculations steps unnecessaryly.

This article has been dead for over six months. Start a new discussion instead.