Changeset 8157

Timestamp:

Mar 20, 2011, 9:42:55 AM (14 years ago)

Author:

steve

Message:

Found annoying bug x(3/2) should be x(1.5)

File:

• trunk/anuga_validation/validation_tests/dam_break_analytical.py (modified) (2 diffs)

trunk/anuga_validation/validation_tests/dam_break_analytical.py

@@ -29 +29 @@
     Taken from formula (3.14) from Zoppou and Roberts SWW notes
     """
-    x = h2/h1
+    
+    h21 = h2/h1
     h01 = h0/h1
     
-    return x**3 -9.0*(x**2)*h01 +16.0*(x**(3/2))*h01 -x*h01*(h01+8.0) +h01**3
+    return h21**3 -9.0*(h21**2)*h01 +16.0*(h21**(1.5))*h01 -h21*h01*(h01+8.0) +h01**3
 
 
 h0 = 1.0
-h1 = 10.0
+h1 = 2.0
+
+h21 = 1.0
+h01 = 0.5
+
+print h21**3 - 9.0*(h21**2)*h01 + 16.0*(h21**(1.5))*h01 - h21*h01*(h01+8.0) + h01**3
+
+h21 = 0.5
+h01 = 0.5
+
+print  - 8.0*h21**3 + 16.0*h21**(2.5) -8*h21**2 
+
 
 n = 10
@@ -54 +66 @@
 print wu(h0), wu(h1)
 
-#from scipy.optimize import brentq
+from scipy.optimize import brentq
 
-#h2 = brentq(wu, h0,h1,disp=True)
+h2 = brentq(wu, h0,h1,disp=True)
 
 #print  h2, wu(h2)