Last change
on this file since 7922 was
7922,
checked in by mungkasi, 13 years ago

Modifying codes so that arrays are compatible with numpy instead of Numeric.

File size:
1.9 KB

Rev  Line  

[7922]  1  from scipy import sin, cos, sqrt, pi, dot 

 2  from gaussPivot import * 

 3  from parameter import * 

 4  

 5  def root_g(a,b,t): 

 6  dx = 0.01 

 7  def g(u): 

 8  return u + 2.0*A*pi/T*sin(2.0*pi/T*(t+u)) 

 9  while 1: 

 10  x1,x2 = rootsearch(g,a,b,dx) 

 11  if x1 != None: 

 12  a = x2 

 13  root = bisect(g,x1,x2,1) 

 14  else: 

 15  break 

 16  return root 

 17  def shore(t): 

 18  a = 0.2#1.0 

 19  b = 0.2#1.0 

 20  #dx = 0.01 

 21  u = root_g(a,b,t) 

 22  xi = 0.5*u*u + A*cos(2.0*pi/T*(t+u)) 

 23  #position = 1.0 + xi 

 24  return [xi, u]#[position, u] 

 25  

 26  

 27  

 28  def prescribe(x,t): 

 29  from Numeric import Float 

 30  from numpy import array, zeros 

 31  if x<1.0: 

 32  q = zeros(2,Float) 

 33  else: 

 34  q = shore(t) 

 35  def fun(q): #Here q=(w, u) 

 36  f = zeros(2,Float) 

 37  f[0] = q[0] + 0.5*q[1]**2.0  A*j0(4.0*pi/T*(1.0+q[0]x)**0.5)*cos(2.0*pi/T*(t+q[1])) 

 38  f[1] = q[1] + A*j1(4.0*pi/T*(1.0+q[0]x)**0.5)*sin(2.0*pi/T*(t+q[1]))/(1+q[0]x)**0.5 

 39  return f 

 40  def newtonRaphson2(f,q,tol=1.0e15): 

 41  for i in range(30): 

 42  h = 1.0e4 ##1.0e4 may be too large. 

 43  n = len(q) 

 44  jac = zeros((n,n),Float) 

 45  if 1.0+q[0]x<0.0: 

 46  #temp1 = 1.0+q[0]x 

 47  #q[0] = q[0]temp1 

 48  q[1] = SomethingWRONG 

 49  #return q 

 50  #q[0] =0.25 

 51  f0 = f(q) 

 52  for i in range(n): 

 53  temp = q[i] 

 54  q[i] = temp + h 

 55  f1 = f(q) 

 56  q[i] = temp 

 57  jac[:,i] = (f1  f0)/h 

 58  if sqrt(dot(f0,f0)/len(q)) < tol: return q 

 59  dq = gaussPivot(jac,f0) 

 60  q = q + dq 

 61  if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q 

 62  print 'Too many iterations' 

 63  q = newtonRaphson2(fun,q) 

 64  return q[0], q[1] 

 65  

 66  #w,u=prescribe(0.0,0.0) 

Note: See
TracBrowser
for help on using the repository browser.