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.

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.