# How to incorporate boundary conditions?

I have this linearized PDE problem:

$$\lambda a^\ast+\bar{a}v^\ast_z+\bar{a}_zv^\ast+\bar{v}a^\ast_z+\bar{v}_za^\ast=0$$

$$\lambda \bar a h^\ast -\bar a \bar h v^\ast_z-\bar a \bar v_z h^\ast - \bar h \bar v_z a^\ast + \bar v \bar h a^\ast_z+ \bar v \bar a_z h^\ast + \bar h \bar a_z v^\ast= -\frac{ F h^\ast}{2}$$

$$\lambda(\bar{h}^2a^\ast-2\bar{a}\bar{h}h^\ast)-2\bar{a}\bar{v}\bar{h}h^\ast_z-2\bar{a}\bar{v}\bar{h}_z h^\ast-2\bar{a}\bar{h}\bar{h}_zv^\ast-2\bar{v}\bar{h}\bar{h}_za^\ast+ \bar{v}\bar{h}^2 a^\ast_z + 2\bar{v}\bar{a}_z \bar{h}h^\ast + \bar{h}^2 \bar{a}_z v^\ast =2{B}\bar{a}a^\ast$$ where $$\bar a$$, $$\bar v$$, and $$\bar h$$ are steady state solutions, $$a^\ast$$, $$v^\ast$$, and $$h^\ast$$ are perturbed quantity, $$F$$ and $$B$$ are parameters and $$\lambda$$ is the complex eigenvalue. The domain for $$z$$ is $$[0,L]$$, here L=10.

The following codes provides eigenvalues of the above problem which has been discretized through the utilization of Cebyshev collocation method. My question is how to change the codes to incorporate boundary conditions: $$a^\ast(0)=0$$, $$v^\ast(0)=0$$ and $$h^\ast(0)=0$$. Please help. Thanks in advance.

ClearAll["Global*"];
ClearAll[recurrenceTables];
recurrenceTables[x_?VectorQ, n_] := Module[{T, Tprime}, T = ConstantArray[0., {n + 1}];
Tprime = ConstantArray[0., {n + 1}];
T[[1]] = ConstantArray[1., Length[x]];
T[[2]] = (2.0 x)/L - 1;
Tprime[[1]] = ConstantArray[0., Length[x]];
Tprime[[2]] = ConstantArray[2.0/L, Length[x]];
Do[T[[i]] = 2.0 ((2.0 x)/L - 1) T[[i - 1]] - T[[i - 2]];
Tprime[[i]] = (4.0/L) T[[i - 1]] +
2.0 ((2.0 x)/L - 1) Tprime[[i - 1]] - Tprime[[i - 2]];, {i, 3,
n + 1}];
{Transpose[T], Transpose[Tprime]}];

n = 300;(*Collocation points*)
(*parameters*)
Br = 1.2; Dr = 2.0; L = 10.0;
B = (Log[Dr] (1.0 - Br^2*Dr))/(L*(1 - Dr) Br^2);
ita = 1.0;
F = (4.0*Log[Dr])/L;

(*Definig Collocation points*)
y = L/2 (Cos[N[Pi/n] Range[0., n]] + 1);
{T, Tprime} = recurrenceTables[y, n];

(*Defining all the functions we need*)
Z[x_] := x;
v[Z_] := Exp[(F/4.0 ita) Z];
a[Z_] := Exp[-(F/4.0 ita) Z];
h[Z_] := Exp[-(F/4.0 ita) Z] Sqrt[((4.0 B)/F) + (1 - ((4.0 B)/F)) Exp[(F/4.0 ita) Z]];

Zy = Map[Z, y];
Vbar = DiagonalMatrix[SparseArray[Map[v, Zy]]];
Vbarprime = DiagonalMatrix[SparseArray[Map[v', Zy]]];
Abar = DiagonalMatrix[SparseArray[Map[a, Zy]]];
Abarprime = DiagonalMatrix[SparseArray[Map[a', Zy]]];
Hbar = DiagonalMatrix[SparseArray[Map[h, Zy]]];
Hbarprime = DiagonalMatrix[SparseArray[Map[h', Zy]]];

(*Defining all the entries for the final U and W matrices*)
T2 = Abar . T; T31 = Hbar . Hbar . T;
T32 = 2.0 Abar . Hbar . T;
U = ArrayFlatten[{{T, 0. T, 0. T}, {0. T, T2, 0. T}, {T31, -T32,
0. T}}];
Tprime = Tprime/L;
A1 = 2.0 Vbar . Tprime + Vbarprime . T;
V1 = 2.0 Abar . Tprime + Abarprime . T;
V2 = Hbar . Abarprime . T - 2.0 Abar . Hbar . Tprime;
A2 = 2.0 Vbar . Hbar . Tprime - Hbar . Vbarprime . T;
H2 = Vbar . Abarprime . T - Abar . Vbarprime . T + (F/(2.0*ita)) T;
H3 = -4.0 Abar . Vbar . Hbar . Tprime - 2.0 Abar . Vbar . Hbarprime . T + 2. Vbar . Abarprime . Hbar . T;
V3 = Hbar . Hbar . Abarprime . T - 2.0 Abar . Hbar . Hbarprime . T;
A3 = 2.0 Vbar . Hbar . Hbar . Tprime - 2.0 Vbar . Hbar . Hbarprime . T - (2.0 B/ita) Abar . T;
W = ArrayFlatten[{{-A1, 0. T, -V1}, {-A2, -H2, -V2}, {-A3, -H3, -V3}}];
(* finding eigenvalues*)
Vals = Union[Eigenvalues[{W, U}]]

• Is it same model as discussed on mathematica.stackexchange.com/questions/289171/… ? Commented Oct 14, 2023 at 6:05
• Yes, it is the same model @AlexTrounev Commented Oct 16, 2023 at 3:33
• There are no a*, v*, h* in your code. How you define these functions? Commented Oct 16, 2023 at 5:32
• a*=TA, v*=TV, h*=TH in matrix form (approximated by chebyshev polynomals)@AlexTrounev Commented Oct 16, 2023 at 19:37
• See my answer. Please pay attention that ChebyshevT[n. x] are orthogonal on $(-1,1)$ only, Check for example Table[2/Pi Integrate[ ChebyshevT[n, x] ChebyshevT[m, x]/Sqrt[1 - x^2], {x, -1, 1}], {n, 3}, {m, 3}]. With this in mind, it is better to map the interval $(0, L)$ to the $(-1,1)$ as in my answer on mathematica.stackexchange.com/questions/289171/… Commented Oct 17, 2023 at 3:38

Here we can use second code from my answer here. Please pay attention that we add boundary conditions to the main system of equations using additional parameters c1,c2,c3.

ClearAll["Global*"];
n = 9;(*Collocation points*)
B = 0.11; ita = 1.0; F = 0.28; L = 10;(*parameters*)

(*Definig Collocation points*)
y = -Cos[N[Pi/n] Range[0., n]];

vec[x_] := Table[ChebyshevT[n, x], {n, 0, 9}]; Map[vec, y];
vec1[x_] := Table[i ChebyshevU[-1 + i, x], {i, 0, 9}]; Map[vec1, y];

(*Defining all the functions we need*)
Z[x_] := L (x + 1)/2;
v[Z_] := Exp[(F/4.0 ita) Z];
a[Z_] := Exp[-(F/4.0 ita) Z];
h[Z_] :=
Exp[-(F/4.0 ita) Z] Sqrt[((4.0 B)/
F) + (1 - ((4.0 B)/F)) Exp[(F/4.0 ita) Z]];

Zy = Map[Z, y];
V0 = DiagonalMatrix[SparseArray[Map[v, Zy]]];
V0p = DiagonalMatrix[SparseArray[Map[v', Zy]]];

A0 = DiagonalMatrix[SparseArray[Map[a, Zy]]];
A0p = DiagonalMatrix[SparseArray[Map[a', Zy]]];

H0 = DiagonalMatrix[SparseArray[Map[h, Zy]]];
H0p = DiagonalMatrix[SparseArray[Map[h', Zy]]];

AA1 = Array[aa, n + 1]; HH1 = Array[hh, n + 1]; VV1 =
Array[vv, n + 1]; a1[y_] := AA1 . vec[y] + c1;
h1[y_] := HH1 . vec[y] + c2;
v1[y_] := VV1 . vec[y] + c3; a1p[y_] := AA1 . vec1[y];
h1p[y_] := HH1 . vec1[y]; v1p[y_] := VV1 . vec1[y];
A1 = Map[a1, y]; H1 = Map[h1, y]; V1 = Map[v1, y]; A1p =
2/L Map[a1p, y]; H1p = 2/L Map[h1p, y]; V1p = 2/L Map[v1p, y];

eq1 = \[Lambda] a1[z] + a[z] v1'[z] + a'[z] v1[z] + v[z] a1'[z] +
v'[z] a1[z];
eq2 = \[Lambda] a[z] h1[z] - a[z] h[z] v1'[z] - a[z] v'[z] h1[z] -
h[z] v'[z] a1[z] + v[z] h[z] a1'[z] + v[z] a'[z] h1[z] +
h[z] a'[z] v1[z] + (F/(2.0*ita)) h1[z];
eq3 = \[Lambda] (h[z]^2 a1[z] - 2.0 a[z] h[z] h1[z]) -
2.0 a[z] v[z] h[z] h1'[z] - 2.0 a[z] v[z] h'[z] h1[z] -
2.0 a[z] h[z] h'[z] v1[z] - 2.0 v[z] h[z] h'[z] a1[z] +
v[z] h[z]^2 a1'[z] + 2.0 v[z] a'[z] h[z] h1[z] +
h[z]^2 a'[z] v1[z] - ((2.0*B)/ita) a[z] a1[z];
bc = {a1[-1] == 0, h1[-1] == 0, v1[-1] == 0};

eqs = Join[\[Lambda] A1 + A0 . V1p + A0p . V1 + V0 . A1p +
V0p . A1, \[Lambda] A0 . H1 - A0 . H0 . V1p - A0 . V0p . H1 -
H0 . V0p . A1 + V0 . H0 . A1p + V0 . A0p . H1 +
H0 . A0p .
V1 + (F/(2.0*ita)) H1, \[Lambda] (H0 . H0 . A1 -
2.0 A0 . H0 . H1) - 2.0 A0 . V0 . H0 . H1p -
2.0 A0 . V0 . H0p . H1 - 2.0 A0 . H0 . H0p . V1 -
2.0 V0 . H0 . H0p . A1 + V0 . H0 . H0 . A1p +
2.0 V0 . A0p . H0 . H1 +
H0 . H0 . A0p . V1 - ((2.0*B)/ita) A0 . A1, bc];

(*creating entries of the final U and W matrix below*)var =
Join[AA1, HH1, VV1, {c1, c2, c3}]; {v, mat} =
CoefficientArrays[eqs, var];

W = Normal[mat] /. \[Lambda] -> 0; U =
W - Normal[mat] /. \[Lambda] -> 1;

{val, fun} = Eigensystem[{W, U}];

val // Union

(* Out[]= {-7.43598 + 0. I, -1.99299 + 0.195764 I, -1.99299 -
0.195764 I, -1.44718 + 0.561469 I, -1.44718 -
0.561469 I, -0.980962 - 1.23891 I, -0.980962 +
1.23891 I, -0.866472 + 0.912756 I, -0.866472 -
0.912756 I, -0.283814 - 0.635496 I, -0.283814 +
0.635496 I, -0.0122723 - 0.0265819 I, -0.0122723 + 0.0265819 I,
0.0799315 - 0.507003 I, 0.0799315 + 0.507003 I, 0.207337 - 2.57139 I,
0.207337 + 2.57139 I, 0.468133 - 0.194187 I,
0.468133 + 0.194187 I, Indeterminate, ComplexInfinity}*)