Javascript required
Skip to content Skip to sidebar Skip to footer
Home / Examples Where Computer Fails to Find Solution Numerical Analysis

Examples Where Computer Fails to Find Solution Numerical Analysis


1.1 Linear equations

Solving linear systems of equations is straightforward using the numpy submodule linalg.solve

Example 1: Solve the following system of linear equations

\left\{\begin{array}{l}x + 3y + 5z = 10\\2x + 5y + z = 8\\2x + 3y + 8z = 3 \end{array} \right.

The exact solution is

\left[\begin{array}{l}x\\ y\\ z\end{array} \right]= \left[ \begin{array}{lll} 1 &3 & 5 \\ 2 & 5 & 1\\ 2& 3& 8 \end{array} \right]^{-1} \left[\begin{array}{l}10\\ 8\\ 3\end{array} \right]= \frac{1}{25} \left[\begin{array}{l}-232\\ 129\\ 19\end{array} \right]= \left[\begin{array}{l}-9.28\\ 5.16\\ 0.76\end{array} \right]

Here is the numerical solution using the function solve of the numpy submodule linalg 

                                    
>>> import numpy as np >>> import np.linalg  >>> A = np.mat('[1,3,5; 2,5,1; 2,3,8]') >>> b = np.mat('[10;8;3]') >>> np.linalg.det(A) -25.000000000000004 The determinant of A is not 0, so there is a unique solution >>> np.linalg.solve(A, b) matrix([[-9.28],         [ 5.16],         [ 0.76]])

1.2 Root finding

To find the roots of a polynomial, used the numpy roots subroutine.

Example 2: Find the roots of the polynomial x^3-3x^2+2x - 1 

                                    >>> import numpy as np >>> p = [1,-3,2,-1] >>> np.roots(p) array([ 2.32471796+0.j, 0.33764102+0.56227951j, 0.33764102-0.56227951j])

To find the roots of a non-linear equations, use the bissection method implemented in the scipy submodule optimize.bisect or the Newton-Raphson method implemented in the scipy submodule optimize.newton .

Example 3. Find the roots of the non-linear equation x + 2cos(x) = 0

  • Bissection method starting on the interval [-2, 2] 
                                                    >>> import scipy >>> import scipy.optimize >>> def f(x):          y = x + 2*scipy.cos(x)          return y           >>> scipy.optimize.bisect(f, -2, 2)        # starting interval [-2, 2] -1.0298665293221347
  • Newton Raphson method with starting point at x0 = 2. 
                                                    >>> import scipy >>> import scipy.optimize >>> def f(x):          y = x + 2*scipy.cos(x)          return y           >>> scipy.optimize.newton(f,  2)      # starting point x0 = 2, xn+1 = xn - f(xn)/f'(xn) for n>0 -1.0298665293222757

To find a root of a set of non-linear equations, the scipy submodule optimize.fsolve is needed.

Example 4: First find the roots of the single-variable non-linear equation x + 2cos(x) = 0 using fsolve at starting point  x0 = 0.3 

                                    >>> import scipy >>> import scipy.optimize  >>> def f(x):          y = x + 2*cos(x)          return y  >>> x0 = scipy.optimize.fsolve(f, 0.3)       # 0.3 is the starting point  >>> print x0, f(x0) -1.0298665293

Example 5: Find the solution of the system of non-linear equations starting at x0 = 1 and x1 = 1 using fsolve

\left\{\begin{array}{l} x_0cos(x_1) = 4\\ x_0x_1 - x_1 = 5 \end{array} \right. 

                                    >>> import scipy >>> import scipy.optimize  >>> def func2(x):         y = [x[0]*cos(x[1]) - 4, x[1]*x[0] - x[1] - 5]         return y  >>> x0 = scipy.optimize.fsolve(func2, [1, 1]) >>> print x0 [ 6.5041,  0.9084]

Python offers an alternative way of defining a function using the lambda form. The lambda form allows to create a function object.

Example 6: Solve the system on non-linear equations 

starting at x=1, y = -1, z =2

\left\{\begin{array}{lll} x^2 -2xy & = & z\\ xy + y^2z & = & 0\\ xz - yz+ x &=& 1 \end{array}\right. 

                                    >>> import scipy >>> import scipy.optimize >>> f = lambda x: [x[0]**2 - 2*x[0]*x[1] - x[2], x[0]*x[1] + x[1]**2*x[2], x[2]*x[0] - x[1]*x[2]+ x[0] - 1]     >>> x0 = scipy.optimize.fsolve(f, [1, -1, 2]) >>> print x0 [ 0.38196601, -0.61803399,  0.61803399]

1.3 Integration

The integrate function provides several integration techniques including an ordinary differential equation integrator.

1.3.1 Integration of Polynomials

Integrate the following polynomial forn 0 to 1 using the numpy function pol.integ
p = (x^5 -4x^3 -2x^2 -3x +2)^5 

                                    >>> import numpy as np >>> pol = np.poly1d([1, 0, -4, -2, -3, 2])**5       # obtain the numpy polynomial form for the integrant >>> polint = pol.integ()                            # close form of the integral >>> print polint(1) - polint(0)                     # apply the integral boundaries to compute the value of the integral -579.48595177

1.3.2 Numerical integration: the scipy.integrate module 

                                    >>> import scipy as sp >>> import sp.integrate >>> help(sp.integrate)               # get help  >>> dir(sp.integrate)                # list the functions available  Methods for Integrating Functions  odeint                               #Integrate ordinary differential equations. quad                                 #General purpose integration. dblquad                              #General purpose double integration. tplquad                              #General purpose triple integration. gauss_quad                           #Integrate func(x) using Gaussian quadrature of order n.

The function scipy.integrate.quad provides a mean to integrate a function of one variable between two points. The points can be \pm \infty  (\pm integrate.inf indicates infinite limits).

Example 7: Compute A = \int_0^2 x^2 \;dx. The exact value of A is: A = \left[ \frac{1}{3}x^3\right]_0^2=\frac{8}{3}\simeq 2.6666.. Here are the steps to find the numerical solution 

                                    >>> import scipy.integrate >>> def func(x):          y = x**2         return y  >>> result, error = scipy.integrate.quad(func, 0.0, 2)  #numerical solution and error >>> result  2.666666666666667 >>> error 2.960594732333751e-14)                                  # the error in negligible

Example 8: Compute A = \int_0^\infty e^-{\frac{1}{2}x}\; dx. The exact solution is 2 

                                    >>> import scipy.integrate >>> import math >>> y = lambda x: math.exp(-0.5*x)                 # using the lambda form to define a (pointer to a) function  >>> result, error = scipy.integrate.quad(y, 0, scipy.inf) >>> result 2.0 >>> error 7.161469913133811e-11

Example 9: Suppose you wish to integrate the bessel function jv(2.5, x) along the interval [0, 4.5].

I = \int_0^{4.5} \; J_{2.5}(x)\; dx 

                                    >>> import scipy.integrate >>> import scipy.special  >>> y = lambda x: scipy.special.jv(2.5,x) >>> result, error = scipy.integrate.quad(y, 0, 4.5) >>> print result 1.1178179380783249 >>> print error 7.8663172481899801e-09)

1.3.3 Double integration 

                                    dblquad                              # double integration

Example 10:

Compute the integral \int_1^{2} \int_0^{1} \frac{x^2}{y} dx dy.   The exact value is \frac{1}{3} log(2)= 0.23104906018664842 

                                    DoubleIntegrationDemo.py """double integration"""  from __future__ import division import scipy import scipy.integrate  def f(x,y): 	s = x**2/y 	return s  if __name__ == '__main__':  	result, error = scipy.integrate.dblquad(lambda y , x: f(x,y),  0, 1, lambda y: 1, lambda y: 2)  	 	print '\n' 	print 'Result =',result  	print 'Error =', error

Result = 0.231049060187
Error = 2.56515986437e-15

Example 10:

Compute the integral \int_0^{\infty} \int_1^{\infty} \frac{e^{-xt}}{t^3}dt dx . The exact value is 1/3 

                                    DoubleIntegration.py """ Double integration """  from __future__ import division import scipy import scipy.integrate  def func(x,t): 	y = scipy.exp(-x*t) / t**3 	return y  if __name__ == '__main__':  	 	result, error = scipy.integrate.dblquad(lambda t, x: func(x,t), 0, scipy.Inf, lambda x: 1, lambda x: scipy.Inf)  	 	print '\n' 	print 'Result =',result  	print 'Error =', error

Result = 0.33333333325
Error = 2.86040699212e-09

1.4 Ordinary differential equations: the scipy.integrate.odeint function

                      odeint                               #Integrate ordinary differential equations.

Example 11: Exponential decay.  Solve the ordinary linear equation \frac{dx}{dt} = -x with initial condition x(0)= 2. The exact solution  is x(t) = 2e^{-t} 
                                          OdeDemo1.py  """ Numerical solution of  ordinary differential equation """  # import modules for solving  import scipy import scipy.integrate  # import module for plotting import pylab as pl  def dx_dt(x, t=0):         # specify the initial time point t0=0 	y = -x 	return y  if __name__ == "__main__":  	t = scipy.linspace(0, 6, 1000)                         # create 1000 time points stating at t0=0 	x0 = 2                                                 # initial condition 	x = scipy.integrate.odeint(dx_dt, x0, t)               # solve using odeint 	 	 	pl.plot(t, x) 	pl.show()

Example 12: Exact versus approximated solution. Solve the non-linear first order equation \frac{dy}{dt} = (1-2t)y^2 with boundary condition y_0 = 1.

Plot the curves of the exact solution f(t) = \frac{1}{t^2-t+1} and the approximated solution obtained with the function odeint on the same graph. 
                                          OdeDemo2.py  """ Numerical solution of  ordinary differential equation. Exact versus approximated solution"""  # import modules for solving  import scipy import scipy.integrate  # import module for plotting import pylab as pl  def dy_dt(y, t=0):         # specify the initial time point t0=0 	z = (1 - 2*t)*y**2 	return z  def f_exact(t):              #exact solution starting at y0=1 	y = 1.0/(t**2 -t+1) 	return y  if __name__ == "__main__":  	t = scipy.linspace(0, 10, 1000)                         # create 1000 time points stating at t0=0 	y0 = 1                                                   # initial condition 	y_approx = scipy.integrate.odeint(dy_dt, y0, t)          # solve using odeint 	f = scipy.vectorize(f_exact) 	y_exact = f(t) 	 #	print len(y_approx), len(y_exact) #	print scipy.shape(y_exact), scipy.shape(y_approx[:,0]) 	 	pl.plot(t, y_approx[:,0], 'ro', t,y_exact , 'bx') 	pl.legend(['Approximated solution','Exact solution']) 	 	pl.savefig('ExactVersusApprox.png') 	pl.show()

Example 13: System of non-linear first order differential equations. Solve the following system non-linear first order Lokta Volterra equations

\left\{\begin{array}{ll} \frac{dx}{dt} & =& x -0.1xy \\ \frac{dy}{dt} & =& -\frac{2}{5}y + 0.7xy \end{array}\right.

with boundary conditions x0 = 10, y0 = 5. Plot the curves of x(t) and y(t) on the same graph for t in the interval [0,15]. 

              
                                                OdeDemo3.py  """ Numerical solution of  system of first order ordinary differential equations """  # import modules for solving  import scipy import scipy.integrate  # import module for plotting import pylab as pl  def dX_dt(X, t=0):         # specify the initial time point t0=0 	y = scipy.array([ X[0] -   0.1*X[0]*X[1] , -2./5*X[1] + 0.7*X[0]*X[1] ]) 	return y  if __name__ == "__main__":  	t = scipy.linspace(0, 15, 1000)                         # create 1000 time points stating at t0=0  	X0 = scipy.array([10, 5])                                    # initials conditions: x0=10  and y0=5  	X, infodict = scipy.integrate.odeint(dX_dt, X0, t, full_output=True) #	infodict['message']                                   #'Integration successful.'  	pl.plot(t, X[:,0],'ro',  t, X[:,1], 'bx')     # x = X[:,0] and y = X[:,1] 	pl.legend(['x(t)','y(t)']) 	pl.plot(X[:,0], X[:,1])           # plot phase portrait 	pl.show()

Examples Where Computer Fails to Find Solution Numerical Analysis

Source: https://www.sites.google.com/a/aims-senegal.org/scipy/roots-finding-numerical-integrations-and-differential-equations