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 = -1.0 |
---|
19 | b = 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)) #This is the displacement |
---|
23 | #position = 1.0 + xi |
---|
24 | return [xi, u]#[position, u] |
---|
25 | |
---|
26 | def w_at_O(t): |
---|
27 | return eps*cos(2.0*pi*t/T) |
---|
28 | |
---|
29 | def u_at_O(t): |
---|
30 | a = -1.01 |
---|
31 | b = 1.01 |
---|
32 | dx = 0.01 |
---|
33 | w = w_at_O(t) |
---|
34 | def fun(u): |
---|
35 | return u + A*j1(4.0*pi/T*(1.0+w)**0.5)*sin(2.0*pi/T*(t+u))/(1.0+w)**0.5 |
---|
36 | while 1: |
---|
37 | x1,x2 = rootsearch(fun,a,b,dx) |
---|
38 | if x1 != None: |
---|
39 | a = x2 |
---|
40 | root = bisect(fun,x1,x2,1) |
---|
41 | else: |
---|
42 | break |
---|
43 | return root |
---|
44 | |
---|
45 | |
---|
46 | def prescribe(x,t): |
---|
47 | from Numeric import Float |
---|
48 | from numpy import array, zeros |
---|
49 | q_shore = shore(t) |
---|
50 | q = zeros(2,Float) |
---|
51 | if x<1.0-abs(q_shore[0]): |
---|
52 | q[0] = w_at_O(t) |
---|
53 | q[1] = u_at_O(t) |
---|
54 | else: |
---|
55 | q = q_shore |
---|
56 | def fun(q): #Here q=(w, u) |
---|
57 | f = zeros(2,Float) |
---|
58 | 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])) |
---|
59 | 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 |
---|
60 | return f |
---|
61 | def newtonRaphson2(f,q,tol=1.0e-15): |
---|
62 | for i in range(30): |
---|
63 | h = 1.0e-4 ##1.0e-4 may be too large. |
---|
64 | n = len(q) |
---|
65 | jac = zeros((n,n),Float) |
---|
66 | if 1.0+q[0]-x<0.0: |
---|
67 | #temp1 = 1.0+q[0]-x |
---|
68 | #q[0] = x-1.0+0.01 |
---|
69 | #print "1.0+q[0]-x=",1.0+q[0]-x |
---|
70 | print "SomethingWRONG" |
---|
71 | #return q |
---|
72 | #q[0] =0.25 |
---|
73 | f0 = f(q) |
---|
74 | for i in range(n): |
---|
75 | temp = q[i] |
---|
76 | q[i] = temp + h |
---|
77 | f1 = f(q) |
---|
78 | q[i] = temp |
---|
79 | jac[:,i] = (f1 - f0)/h |
---|
80 | if sqrt(dot(f0,f0)/len(q)) < tol: return q |
---|
81 | dq = gaussPivot(jac,-f0) |
---|
82 | q = q + dq |
---|
83 | if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q |
---|
84 | print 'Too many iterations' |
---|
85 | q = newtonRaphson2(fun,q) |
---|
86 | return q[0], q[1] |
---|
87 | |
---|
88 | #w,u=prescribe(0.0,0.0) |
---|