1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
|
ex8.py: Nonlinear eigenproblem with split form
==============================================
This example solves a nonlinear eigenvalue problem where the
nonlinear function is expressed in split form.
We want to solve the following parabolic partial differential
equation with time delay :math:`\tau`
.. math::
u_t &= u_{xx} + a u(t) + b u(t-\tau) \\
u(0,t) &= u(\pi,t) = 0
with :math:`a = 20` and :math:`b(x) = -4.1+x (1-e^{x-\pi})`.
Discretization leads to a DDE of dimension :math:`n`
.. math::
-u' = A u(t) + B u(t-\tau)
which results in the nonlinear eigenproblem
.. math::
(-\lambda I + A + e^{-\tau\lambda}B)u = 0.
The full source code for this demo can be `downloaded here
<../_static/ex8.py>`__.
Initialization is similar to previous examples. In this case we also
need to import some math symbols.
::
import sys, slepc4py
slepc4py.init(sys.argv)
from petsc4py import PETSc
from slepc4py import SLEPc
from numpy import exp
from math import pi
Print = PETSc.Sys.Print
This script has two command-line options: the discretization size ``n``
and the time delay ``tau``.
::
opts = PETSc.Options()
n = opts.getInt('n', 128)
tau = opts.getReal('tau', 0.001)
a = 20
h = pi/(n+1)
Next we have to set up the solver. In this case, we are going to
represent the nonlinear problem in split form, i.e., as a sum of
terms made of a constant matrix multiplied by a scalar nonlinear
function.
::
nep = SLEPc.NEP().create()
The first term involves the identity matrix.
::
Id = PETSc.Mat().createConstantDiagonal([n, n], 1.0)
The second term has a tridiagonal matrix obtained from the
discretization, :math:`A = \frac{1}{h^2}\operatorname{tridiag}(1,-2,1) + a I`.
::
A = PETSc.Mat().create()
A.setSizes([n, n])
A.setFromOptions()
rstart, rend = A.getOwnershipRange()
vd = -2.0/(h*h)+a
vo = 1.0/(h*h)
if rstart == 0:
A[0, :2] = [vd, vo]
rstart += 1
if rend == n:
A[n-1, -2:] = [vo, vd]
rend -= 1
for i in range(rstart, rend):
A[i, i-1:i+2] = [vo, vd, vo]
A.assemble()
The third term includes a diagonal matrix :math:`B = \operatorname{diag}(b(x_i))`.
::
B = PETSc.Mat().create()
B.setSizes([n, n])
B.setFromOptions()
rstart, rend = B.getOwnershipRange()
for i in range(rstart, rend):
xi = (i+1)*h
B[i, i] = -4.1+xi*(1.0-exp(xi-pi));
B.assemble()
B.setOption(PETSc.Mat.Option.HERMITIAN, True)
Apart from the matrices, we have to create the functions, represented with
`FN` objects: :math:`f_1=-\lambda, f_2=1, f_3=\exp(-\tau\lambda)`.
::
f1 = SLEPc.FN().create()
f1.setType(SLEPc.FN.Type.RATIONAL)
f1.setRationalNumerator([-1, 0])
f2 = SLEPc.FN().create()
f2.setType(SLEPc.FN.Type.RATIONAL)
f2.setRationalNumerator([1])
f3 = SLEPc.FN().create()
f3.setType(SLEPc.FN.Type.EXP)
f3.setScale(-tau)
Put all the information together to define the split operator. Note that
``A`` is passed first so that `SUBSET <petsc4py.PETSc.Mat.Structure.SUBSET>`
nonzero_pattern can be used.
::
nep.setSplitOperator([A, Id, B], [f2, f1, f3], PETSc.Mat.Structure.SUBSET)
Now we can set some options and call the solver.
::
nep.setTolerances(tol=1e-9)
nep.setDimensions(1)
nep.setFromOptions()
nep.solve()
Once the solver has finished, we print some information together with
the computed solution. For each computed eigenpair, we print the
eigenvalue and the residual norm.
::
its = nep.getIterationNumber()
Print("Number of iterations of the method: %i" % its)
sol_type = nep.getType()
Print("Solution method: %s" % sol_type)
nev, ncv, mpd = nep.getDimensions()
Print("")
Print("Subspace dimension: %i" % ncv)
tol, maxit = nep.getTolerances()
Print("Stopping condition: tol=%.4g" % tol)
Print("")
nconv = nep.getConverged()
Print( "Number of converged eigenpairs %d" % nconv )
if nconv > 0:
x = Id.createVecs('right')
x.set(1.0)
Print()
Print(" k ||T(k)x||")
Print("----------------- ------------------")
for i in range(nconv):
k = nep.getEigenpair(i, x)
res = nep.computeError(i)
if k.imag != 0.0:
Print( " %9f%+9f j %12g" % (k.real, k.imag, res) )
else:
Print( " %12f %12g" % (k.real, res) )
Print()
|
