Last change
on this file since 7922 was
7922,
checked in by mungkasi, 15 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.0e-15): |
---|
| 41 | for i in range(30): |
---|
| 42 | h = 1.0e-4 ##1.0e-4 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.