Labels

Showing posts with label PYTHON. Show all posts
Showing posts with label PYTHON. Show all posts

Friday, December 11, 2020

A simple guide to Install scikit ODES

Here is a simple guid to install scikit-odes. Odes is a scikit toolkit for scipy to add extra ode solvers. Specifically it interfaces the Sundials solvers cvode, cvodes, ida and idas.


a) Install BLAS and LAPACK

>sudo apt-get install libopenblas-dev liblapack-dev


b) Install Sundials 5.1.0

>wget https://computing.llnl.gov/projects/sundials/download/sundials-5.1.0.tar.gz

>tar xvf sundials-5.1.0.tar.gz 

>mkdir builder 

>cd builder 

>cmake ../../sundials-5.1.0/ -DLAPACK_ENABLE=ON 

>make 

>sudo make install

c) Install scikits odes

add export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/bin to .bashrc 

>source ~/.bashrc 

>sudo pip install scikits.odes


Reference:

https://scikits-odes.readthedocs.io/en/latest/installation.html

Thursday, August 27, 2020

Scipy optimize error "length of x0 != length of bounds"

 

from scipy.optimize import minimize
def fun(x):
    return x-1
minimize(fun, [2],bounds = ((0,4)))
It will raise "ValueError: length of x0 != length of bounds". We should change ((0,4)) to ((0,4),)
minimize(fun, [2], bounds=((0, 4),))
Reference: https://github.com/scipy/scipy/issues/9692#issuecomment-455258618

Tuesday, August 25, 2020

Compute Jacobian matrix by Sympy

Computing the Jacobian matrix is not easy, but we can use the python symbolical package Sympy to get it.

Here is a simple example.

import sympy as sym
from sysmpy import init_printing
from sympy import sin,cos,Matrix

init_printing()

x,y=sym.symbols('x y')
F=Matrix([x*cos(y),x*sin(y),x**2])
X=Matrix([x,y])
print(F.jacobian(X))
print(F.jacbian(X).subs([(x,0),(y,1)]))



Reference:
https://docs.sympy.org/latest/modules/matrices/matrices.html
https://www.sympy.org/scipy-2017-codegen-tutorial/notebooks/20-ordinary-differential-equations.html

Thursday, June 18, 2020

Runge-Kutta Method

The fourth order Runge-Kutta method (RK4) is the most widely used algorithm for solving an ordinary differential equation.
\[ dy/dx = f(x,y) \]
\[ y_{n+1} = y_n +\frac{h}{6} (k_1+k_2+k_3+k_4) \],
where
\[ k_1 = f(x_n,y_n) \]
\[ k_2 = f(x_n+\frac{h}{2}, y_n+\frac{h}{2} k_1 \]
\[ k_3 = f(x_n+\frac{h}{2}, y_n+\frac{h}{2} k_2 \]
\[ k_4 = f(x_n+h,y_n+hk_3) .\]

$k_1$ is the slope at the beginning of the interval $[x_n,x_n+h]$.
$k_2$ is the slope at the midpoint of the interval, using slope $k_1$ to determine the value of $y$ at the point $x_n+h/2$ using Euler's method. 
$k_3$ is again the slope at the midpoint, but now using the slope $k_2$ to determine the $y$ value, $y(x_0+n*h+h/2) = y(x_0+n*h)+h/2*f(x_0+n*h+h/2,y(x_0+n*h)+h/2*k1)$.
$k_4$ is the slope at the end of the interval, with its $y$ value determined using $k_3$.
  

Tuesday, June 9, 2020

Python notes


import numpy as np
np.power(-1.2,5./3.)
>>>RuntimeWarning: invalid value encountered in power
nan

pow(-1.2,5./3.)
>>>0.6775459407943406-1.173543993917852j
np.seterr(invalid='raise')
try:
	np.power(-1.2,5./3.)
except:
    print('error here')  

def replaceZeroes(data):
  min_nonzero = np.min(data[np.nonzero(data)])
  data[data == 0] = min_nonzero
  return data
  

install numpy for python2.7
sudo apt-get install python-pip # install pip for python2 , will get pip2
/usr/bin/python2.7 -m pip install numpy==1.14.6
python setup.py install --prefix=xxx
Python.h not found
sudo apt-get install python-dev # for python2
sudo apt-get install pythonx-dev # x is version of python, i.e. 3.7
sudo apt-get install python-matplotlib # for python2
sudo apt-get install python-scipy
from scipy.misc import derivative
derivative(exp(x), 1.0, dx=1e-6)
To modify the global variable inside the function, we need add the keyword 'global'
a=0
def fun():
	global a
    a+=1
time
import time
t0=time.time()
time.time() - t0
numpy insert

array=np.insert(array,0,'s')

if any(x>5 for x in range(-2,20)):
	print('one element >5')
fill matrix with same element
np.full((2,2),1)
construct matrix
from scipy.sparse import diags
diagonals = [[1, 2, 3, 4], [1, 2, 3], [1, 2]]
diags(diagonals, [0, -1, 2]).toarray()
diags([1,1],[-1,1],shape=(3,4)).toarray()
print numpy narray nice
def matprint(mat, fmt="g"):
    col_maxes = [max([len(("{:"+fmt+"}").format(x)) for x in col]) for col in mat.T]
    for x in mat:
        for i, y in enumerate(x):
            print(("{:"+str(col_maxes[i])+fmt+"}").format(y), end="  ")
        print("")
print(np.array2string(A, suppress_small=True, formatter={'float': '{:0.4f}'.format}))
pip update package
pip install package-name -U --user


a=[1,2,3,4]
a[1:-1] will give 2,3 doesn't include last element 3!
a[1:4] or a[1:] will give 2,3,4 
a[-1] give 4


a=[2]*number # initialize list with same value 2
b=np.full((1,3),2) #give array([[2, 2, 2]])
b[0] give array([2, 2, 2])
c=np.full((3),2) give array([2, 2, 2])
c[0] give 2

insert one row into numpy array
a=np.zeros((2,2))
#line 1 is the line before to insert, 0 is the axis 
a=np.insert(1,np.array((2,2)),0)
give
array([[ 0.,  0.],
       [ 2.,  2.],
       [ 0.,  0.]])
row_number,col_number = a_array.shape
row_number = len(a_array)
sympy.latex() converts mathematical expressions into Latex equations.

check varible type is numpy.float64
isinstance(numpy.float64(1.3), numpy.float64) 
>>> True
plot odd numbers subplot in matplotlib
fig, axes = plt.subplots(2, 3, figsize=(20, 10))
fig.delaxes(ax=axes(2,2))
we will get 5 subplots. matplotlib subplot tight_layout
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
save arrays to file
import numpy
list1 = [1, 2, 3, 4]
list2 = [0.45, 0.98, 0.89, 0.21]
dat = numpy.array([list1, list2])
dat = dat.T
numpy.savetxt('data.txt', dat, delimiter = ',')



Reference:
https://stackoverflow.com/questions/3419082/write-multiple-numpy-arrays-to-file
https://docs.scipy.org/doc/scipy/reference/generated/scipy.misc.derivative.html
https://gist.github.com/braingineer/d801735dac07ff3ac4d746e1f218ab75
https://stackoverflow.com/questions/8298797/inserting-a-row-at-a-specific-location-in-a-2d-array-in-numpy

Monday, June 8, 2020

Test two float number if equal

math.isclose(a,b,rel_tol=1e-9,abs_tol=0)
rel_tol is the relative tolerance – it is the maximum allowed difference between a and b, relative to the larger absolute value of a or b.
abs_tol is the minimum absolute tolerance
the result will be: abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
Return True if the values a and b are close to each other and False otherwise.

  
  import math
  import sys
  print(math.isclose(0.3,0.3+1e-7,rel_tol=sys.float_info.epsilon,abs_tol=0))

bool CompareDoubles2 (double A, double B) 
{
   diff = A - B;
   return (diff < EPSILON) && (-diff < EPSILON);
}
Reference: https://stackoverflow.com/questions/17333/what-is-the-most-effective-way-for-float-and-double-comparison

Sunday, June 7, 2020

Install scipy 1.4.1

Here is the method to install Scipy 1.4.1.

python3 -m pip install --upgrade pip
pip3 install pybind11
pip3 install --no-cache --upgrade scipy  or pip3 install --upgrade scipy==1.4.1

Reference:

Saturday, June 6, 2020

Numerical solve the ordinary differential equation with python

Some excellent instructions are listed below.
It is real very easy to get it.


 Just do it!

multiple shooting
A First Course in Ordinary Differential Equations (book)

newton method

henyey relaxation

Finite Difference Schemes from 1st to 4th order

Implicit Euler method

SUite of Nonlinear and DIfferential/ALgebraic Equation Solvers

Solving differential equations using modified Picard iteration

sympy

Wednesday, October 10, 2018

Python symbol computation package sympy


SymPy is a Python library for symbolic mathematics. It aims to become a full-featured computer algebra system (CAS) while keeping the code as simple as possible in order to be comprehensible and easily extensible. SymPy is written entirely in Python.

online sympy https://www.sympygamma.com/

from sympy import *

#define variable
x = Symbol('x')
y = Symbol('y')
a, b, c = symbols('a,b,c')

#expand
expand((x + y)**3)
expand(sin(x+y),trig=True)

#simplify
simplify((x + x*y) / x)

#limit
#limit(function, variable, point)
limit(sin(x)/x,x,0)
limit(1/x, x, oo)
limit( (sin(x+y) - sin(x))/y,y,0)

#diff
#diff(func,var,n)
diff(sin(x),x)
diff(sin(x),x,2)

#taylor series
#series(expr, var)
series(sin(x),x)

#integrate
integrate(x**2,x)
integrate(2*x + sinh(x), x)
integrate(exp(-x**2)*erf(x), x)
-----
integrate(x**3, (x, -1, 1))
integrate(exp(-x), (x, 0, oo))
integrate(exp(-x**2), (x, -oo, oo))

#factor
factor(x**2-1)

#solve equation
solve(x**4 - 1, x)
solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])

#matrix
A = Matrix([[1,x], [y,1]])

#dsolve
f, g = symbols('f g', cls=Function)
dsolve(f(x).diff(x, x) + f(x), f(x))

# real and positive variable
from sympy import symbols integrate
a,b,k=symbols('a,b,k',real=True,positive=True)
f=k**2/(k**2+a**2)/(1+k**2*b**2)**(5/6)
g=integrate(f,(k,-oo,oo))


REFERENCE
https://wizardforcel.gitbooks.io/scipy-lecture-notes/content/15.html