|
| 1 | +# Tutorial: Using NEP-PACK from python |
| 2 | + |
| 3 | +## PyJulia |
| 4 | + |
| 5 | +The previous tutorial illustrated how a NEP defined |
| 6 | +in python code can be solved with NEP-PACK. |
| 7 | +Python is currently a more mature language than Julia, |
| 8 | +and there are considerable packages and features |
| 9 | +in python not present in Julia. If you need these features, |
| 10 | +it may be more convenient to call NEP-PACK |
| 11 | +from python, rather than calling python code from julia. |
| 12 | + |
| 13 | +The python package [PyJulia](https://github.com/JuliaPy/pyjulia) |
| 14 | +gives us that possibility. The installation of PyJulia on Ubuntu linux is simple: |
| 15 | +``` |
| 16 | +$ python3 -m pip install julia # Only necessary first time you run it |
| 17 | +... |
| 18 | +$ python3 |
| 19 | +>>> import julia |
| 20 | +>>> jl = Julia(compiled_modules=False) # compilation flag necessary on ubuntu |
| 21 | +>>> julia.install() # Only necessary first time you run it |
| 22 | +>>> from julia import Base |
| 23 | +>>> Base.MathConstants.golden # Julia's definition of golden ratio |
| 24 | +1.618033988749895 |
| 25 | +``` |
| 26 | + |
| 27 | + |
| 28 | +## Using PyJulia and NEP-PACK |
| 29 | + |
| 30 | +The [`Mder_NEP`](@ref)-function provides a convenient |
| 31 | +way to define NEPs by only |
| 32 | +using a function that computes the matrix ``M(λ)`` |
| 33 | +and its derivatives. |
| 34 | +Let us first define a function which does that in python. We consider |
| 35 | +the problem |
| 36 | +```math |
| 37 | +M(λ)=\begin{bmatrix}3&2\newline3&-1\end{bmatrix}+ |
| 38 | +λ\begin{bmatrix}0&2\newline0&1\end{bmatrix}+ |
| 39 | +e^{0.5 λ}\begin{bmatrix}1&1\newline1&1\end{bmatrix} |
| 40 | +``` |
| 41 | +and implement it with this python code: |
| 42 | +```python |
| 43 | +import numpy as np; |
| 44 | +import cmath as m; |
| 45 | +def my_compute_M(s,der): |
| 46 | + """Compute the matrix M^{(k)}(s) for a given eigenvalue approximation and derivative k""" |
| 47 | + A=np.matrix('3 2; 3 -1'); B=np.matrix('0 2; 0 1'); C=np.matrix('1 1; 1 1'); |
| 48 | + tau=0.5; |
| 49 | + M=pow(tau,der)*m.exp(tau*s)*C |
| 50 | + if (der==0): |
| 51 | + M=M+A+s*B; |
| 52 | + elif (der==1): |
| 53 | + M=M+B; |
| 54 | + return M |
| 55 | +``` |
| 56 | +An evaluation of the matrix function can be done by the call: |
| 57 | +``` |
| 58 | +>>> my_compute_M(0.3,0) |
| 59 | +matrix([[ 6.76183424+0.j, 3.16183424+0.j], |
| 60 | + [ 4.16183424+0.j, -3.7 +0.j]]) |
| 61 | +``` |
| 62 | +We instantiate a new NEP based with `Mder_NEP` which first must be imported |
| 63 | +``` |
| 64 | +>>> from julia.NonlinearEigenproblems import Mder_NEP |
| 65 | +>>> n=2; # Size of the problem |
| 66 | +>>> nep=Mder_NEP(2,my_compute_M); |
| 67 | +``` |
| 68 | +and we can apply most of our solvers to this problem by first importing the corresponding function, in this case we use [`contour_beyn`](@ref). |
| 69 | +``` |
| 70 | +>>> from julia.NonlinearEigenproblems import contour_beyn; |
| 71 | +>>> sol=contour_beyn(nep,logger=1,neigs=1,radius=3) |
| 72 | +Computing integrals |
| 73 | +NonlinearEigenproblems.NEPSolver.MatrixTrapezoidal: computing G... |
| 74 | +NonlinearEigenproblems.NEPSolver.MatrixTrapezoidal: summing terms........................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................ |
| 75 | +Computing SVD prepare for eigenvalue extraction p=1 |
| 76 | +Computing eigenvalues |
| 77 | +Computing eigenvectors |
| 78 | +>>> |
| 79 | +``` |
| 80 | +We can verify that we computed a solution as follows |
| 81 | +``` |
| 82 | +>>> s=sol[0][0]; v=sol[1] |
| 83 | +>>> my_compute_M(s,0)*v |
| 84 | +matrix([[1.71634841e-17-1.59872116e-14j], |
| 85 | + [9.55210099e-17-3.99680289e-15j]]) |
| 86 | +>>> from numpy.linalg import norm |
| 87 | +>>> norm(my_compute_M(s,0)*v) |
| 88 | +1.6479526251408437e-14 |
| 89 | +``` |
| 90 | +Note that in order to obtain better efficiency for |
| 91 | +large-scale problems, and reduce overhead, |
| 92 | +you may want to use [`Mder_Mlincomb_NEP`](@ref), |
| 93 | +as described in the [previous tutorial](tutorial_call_python.md). |
| 94 | + |
| 95 | + |
| 96 | + |
0 commit comments