Added a few examples of scipy solvers

Author: Jonathan Rocher

scipy/fsolve.py

@@ -0,0 +1,22 @@
+""" Testing fsolve
+"""
+from numpy import cos, exp, sin, pi
+from scipy.optimize import fsolve
+
+
+def equations(x,a,b,c):
+    x0, x1, x2 = x
+    eqs = \
+    [3 * x0 - cos(x1*x2) + a,
+    x0**2 - 81*(x1+0.1)**2 + sin(x2) + b,
+    exp(-x0*x1) + 20*x2 + c]
+    return eqs
+
+a = -0.5
+b = 1.06
+c = (10 * pi - 3.0) / 3
+# Optimization start location.
+initial_guess = [0.1, 0.1, -0.1]
+# Solve the system of non-linear equations.
+root = fsolve(equations, initial_guess, args=(a, b, c))
+print equations(root, a, b, c)

scipy/linalg_solve.py

@@ -0,0 +1,34 @@
+from numpy import *
+from scipy.linalg import inv, solve, lu_factor, lu_solve, det, lstsq
+from numpy.linalg import cond
+
+A = array([[ 1.        ,  2.67499884,  7.15561882 ],
+       [ 1.        ,  2.44819414,  5.99365454],
+       [ 1.        ,  2.30846879,  5.32902817]], dtype = float64)
+       
+print "Near singular matrix if cond is large or det is close to 0", cond(A), det(A)
+
+inverse_matrix = inv(A)
+if not allclose(A.dot(inverse_matrix), identity(3)):
+    raise ValueError("Matrix inversion in viscosity computation failed.")
+
+Y = array([1e15, 1e3, 1.], dtype = float64)
+#Y = array([5, 3, 1.], dtype = float64)
+
+# Manual way
+X = dot(inverse_matrix,Y)
+
+# More efficient way: the solver doesn't compute the full inverse, since it is not needed.
+X = solve(A,Y)
+
+# LU decomposition actually returns just like solve: it probably implements LU behind the scene.
+lu, piv = lu_factor(A)
+X = lu_solve((lu, piv), Y)
+
+# The least square solver offers a different algorithm to compute this in order
+# to minimize |Y-AX|.
+X, residues, rank, singular_values = lstsq(A, Y)
+
+if not allclose(dot(A,X), Y):
+    raise ValueError("The solution from the inversion is incorrect: %s"
+                      % dot(A,X))

scipy/optimize.py

@@ -0,0 +1,25 @@
+""" Demo scipy curve_fit
+"""
+from scipy.optimize import curve_fit
+from scipy.stats import norm
+
+from numpy import exp, sin, pi, linspace, abs
+
+def function(x, a , b, f, phi):
+    result = a * exp(-b * sin(f * x + phi))
+    return result
+
+# Create a noisy data set.
+actual_params = [3, 2, 1, pi/4]
+x = linspace(0,2*pi,25)
+exact = function(x, *actual_params)
+noisy = exact + 0.3*norm.rvs(size=len(x)) 
+# Use curve_fit to estimate the function parameters from the noisy data.
+initial_guess = [1,1,1,1]
+estimated_params, err_est = curve_fit(function, x, noisy, p0=initial_guess) 
+print "params", estimated_params
+#array([3.1705, 1.9501, 1.0206, 0.7034])
+# err_est is an estimate of the covariance matrix of the estimates
+#  (i.e. how good of a fit is it)
+print "Estimated covariance diagonal:", err_est.diagonal()
+print "Real errors:", abs(actual_params-estimated_params)