I am working on research article with link https://www.sciencedirect.com/science/article/abs/pii/S0378475421004560
It contains dimensionless non-linear equations mentioned in article.
I am writing the dimensionless linear equations after quasilinearization process that are:
$F_{1\eta\eta}+A_1^jF_{1\eta}+A_2^jF_1+A_3^jF_{1\xi}+A_4^jF_2+A_5^jG+A_6^jH+A_7^jS+A_8^jP=U_1^j$
$F_{2\eta\eta}+B_1^jF_{2\eta}+B_2^jF_2+B_3^jF_{2\xi}+B_4^jF_1=U_2^j$
$G_{\eta\eta}+C_1^jG_\eta+C_2^jG_\xi+C_3^jF_1+C_4^jP_\eta=U_3^j$
$H_{\eta\eta}+D_1^jH_\eta+D_2^jH_\xi+D_3^jF_1=U_4^j$
$S_{\eta\eta}+E_1^jS_\eta+E_2^jS_\xi+E_3^jF_1=U_5^j$
$P_{\eta\eta}+I_1^jP_\eta+I_2^jP_\xi+I_3^jF_1+I_4^jG_{\eta\eta}=U_6^j$
I created the difference equations using backward difference in $\xi-$direction and central difference in $\eta-$direction that are:
$\bigg(1-\dfrac{A1\Delta\eta}{2}\bigg)F_1(i,j-1)+(-2+A2\Delta\eta^2)F_1(i,j)+\bigg(1+\dfrac{A1\Delta\eta}{2}\bigg)F_1(i,j+1)+A4\Delta\eta^2F_2(i,j)=\Delta\eta^2(U1-A3F_{1\xi}(i,j)-A5G(i,j)-A6H(i,j)-A7S(i,j)-A8P(i,j))$
$B4\Delta\eta^2F1(i,j)+\bigg(1-\dfrac{B1\Delta\eta}{2}\bigg)F2(i,j-1)+(-2+B2\Delta\eta^2)F2(i,j)+\bigg(1+\dfrac{B1\Delta\eta}{2}\bigg)F2(i,j+1)=\Delta\eta^2(U2-B3F_{2\xi})$
$\bigg(1-\dfrac{C1\Delta\eta^2}{2}\bigg)G(i,j-1)-2G(i,j)+\bigg(1+\dfrac{C1\Delta\eta^2}{2}\bigg)G(i,j+1)=\Delta\eta^2(U3-C2G_\xi(i,j)-C3F_1(i,j)-C4P_\eta(i,j))$
$\bigg(1-\dfrac{D1\Delta\eta^2}{2}\bigg)H(i,j-1)-2H(i,j)+\bigg(1+\dfrac{D1\Delta\eta^2}{2}\bigg)H(i,j+1)=\Delta\eta^2(U4-D2H_\xi(i,j)-D3F_1(i,j))$
$\bigg(1-\dfrac{E1\Delta\eta^2}{2}\bigg)S(i,j-1)-2S(i,j)+\bigg(1+\dfrac{E1\Delta\eta^2}{2}\bigg)S(i,j+1)=\Delta\eta^2(U5-E2H_\xi(i,j)-E3F_1(i,j))$
$\bigg(1-\dfrac{I1\Delta\eta^2}{2}\bigg)P(i,j-1)-2P(i,j)+\bigg(1+\dfrac{I1\Delta\eta^2}{2}\bigg)P(i,j+1)=\Delta\eta^2(U6-I2P_\xi(i,j)-I3F_1(i,j)-I4G_{\eta\eta}(i,j))$
The relation between $F1$ and $f$ is that $F1=\dfrac{\partial f}{\partial \eta}$ and $f_\xi=\dfrac{\partial f}{\partial \xi}$
Actually I am new to Mathematica but I tried on my own to create code for solution of the system of equations mentioned in article but actually now I am facing trouble in $\xi-$stepping of the code. I can't create the tridiagonal matrices very well, I think so.
The code I have created for their solution is:
ClearAll["Global`*"];
(*Varga's Algorithm for block system with error tracking*)
ClearAll[VargaSolve];
VargaSolve[AF1_, AF2_, AG_, AH_, AS_, AP_, B12_, B21_, b1_, b2_, b3_,
b4_, b5_, b6_, F10_, F20_, G0_, H0_, S0_, P0_,(*initial guesses*)
tol_ : 1.*10^-4, maxIter_ : 100, verbose_ : True ] :=
Module[{F1 = F10, F2 = F20, G = G0, H = H0, S = S0, P = P0,
err = \[Infinity], k = 0, resU, resV, resW, resX, resY, resZ, errU,
errV, errW, errX, errY, errZ, errorList = {}, solveF1, solveF2,
solveG, solveH, solveS,
solveP},(*Precompute linear solvers if matrices are constant*)
solveF1 = LinearSolve[AF1];
solveF2 = LinearSolve[AF2];
solveG = LinearSolve[AG];
solveH = LinearSolve[AH];
solveS = LinearSolve[AS];
solveP = LinearSolve[AP];
If[verbose, Print["Iter\tError"]];
While[k < maxIter && err > tol,(*Residuals*)
resU = b1 - (AF1 . F1 + B12 . F2);
resV = b2 - (B21 . F1 + AF2 . F2);
resW = b3 - AG . G;
resX = b4 - AH . H;
resY = b5 - AS . S;
resZ = b6 - AP . P;
(*Update solutions*)
F1 += solveF1[resU];
F2 += solveF2[resV];
G += solveG[resW];
H += solveH[resX];
S += solveS[resY];
P += solveP[resZ];
(*Compute error*)
errU = Norm[resU, \[Infinity]];
errV = Norm[resV, \[Infinity]];
errW = Norm[resW, \[Infinity]];
errX = Norm[resX, \[Infinity]];
errY = Norm[resY, \[Infinity]];
errZ = Norm[resZ, \[Infinity]];
err = Max[errU, errV, errW, errX, errY, errZ];
AppendTo[errorList, err];
If[verbose, Print[k, "\t", ScientificForm[err, 3]]];
k++];
{F1, F2, G, H, S, P, errorList}]
(*Parameters*)
\[CapitalLambda] = 5; (*Casson fluid parameter*)
Pr = 0.7; (*Prandtl number*)
Sc1 = 160; (*Schmidt number for liquid hydrogen*)
Sc2 = 360; (*Schmidt number for liquid ammonia*)
Le = 10; (*Lewis number*)
Nb = 0.1; (*Brownian diffusion parameter*)
Nt = 0.1; (*Thermophoresis parameter*)
Nr = 0.1; (*Nanoparticles buoyancy ratio parameter*)
\[Lambda] = -2; (*Mixed convection parameter*)
\[Chi] = 5; (*Rotational parameter*)
Nc1 = 0.1; (*Buoyancy ratio parameters for liquid hydrogen*)
Nc2 = 0.1; (*Buoyancy ratio parameters for liquid ammonia*)
M = 0.5; (*Magnetic Parameter*)
(*Defining Grid*)
Subscript[\[Xi], 0] = 0; Subscript[\[Xi], max] = 0.5; \
Subscript[\[Eta], 0] = 0; Subscript[\[Eta], max] = 5; N\[Xi] = 4; N\
\[Eta] = 4; d\[Xi] = (Subscript[\[Xi], max] - Subscript[\[Xi], 0])/
N\[Xi]; d\[Eta] = (Subscript[\[Eta], max] - Subscript[\[Eta], 0])/
N\[Eta];
\[Xi]Grid =
Range[Subscript[\[Xi], 0], Subscript[\[Xi], max], d\[Xi]];
\[Eta]Grid =
Range[Subscript[\[Eta], 0], Subscript[\[Eta], max], d\[Eta]];
L1 = Table[1 - Cos[\[Xi]Grid[[i]]], {i, 1, Length[\[Xi]Grid]}];
L2 = Table[1 + Cos[\[Xi]Grid[[i]]], {i, 1, Length[\[Xi]Grid]}];
L3 = Table[2 + Cos[\[Xi]Grid[[i]]], {i, 1, Length[\[Xi]Grid]}];
Print["L1=", L1[[2 ;; -2]]];
Print["L2=", L2[[2 ;; -2]]];
Print["L3=", L3[[2 ;; -2]]];
B = Table[(1/3)*L3[[i]]*(L2[[i]])^-1*Tan[\[Xi]Grid[[i]]/2], {i, 1,
Length[\[Xi]Grid]}];
\[Beta]5 = Table[(1/2)*(L1[[i]])^2*L3[[i]], {i, 1, Length[\[Xi]Grid]}];
\[Beta] =
Table[(2/3)*L3[[i]]*(L2[[i]])^2*Cos[\[Xi]Grid[[i]]], {i, 1,
Length[\[Xi]Grid]}];
s = Table[(8/27)*L3[[i]]*L2[[i]]^(-2), {i, 1, Length[\[Xi]Grid]}];
\[Beta]4 =
Table[L1[[i]]*L2[[i]]^(-1)*L3[[i]], {i, 1, Length[\[Xi]Grid]}];
\[Alpha] = Table[\[Chi]*\[Beta][[i]], {i, 1, Length[\[Xi]Grid]}];
Print["B=", B[[2 ;; -2]]]
Print["\[Beta]5=", \[Beta]5[[2 ;; -2]]]
Print["\[Beta]=", \[Beta][[2 ;; -2]]]
Print["s=", s[[2 ;; -2]]]
Print["\[Beta]4=", \[Beta]4[[2 ;; -2]]]
Print["\[Alpha]=", \[Alpha][[2 ;; -2]]]
(*Initialize solution arrays*)
F1 = ConstantArray[
0.0, {N\[Xi] + 1,
N\[Eta] + 1}]; F2 = F1; G = F1; H = F1; S = F1; P = F1; f = F1;
(*Initial conditions at \[Xi]=0*)
f0[\[Xi]_, \[Eta]_] := \[Eta] - 1 + Exp[-\[Eta]];
F10[\[Xi]_, \[Eta]_] := 1 - Exp[-\[Eta]];
F20[\[Xi]_, \[Eta]_] := Exp[-\[Eta]];
G0[\[Xi]_, \[Eta]_] := Exp[-\[Eta]];
H0[\[Xi]_, \[Eta]_] := Exp[-\[Eta]];
S0[\[Xi]_, \[Eta]_] := Exp[-\[Eta]];
P0[\[Xi]_, \[Eta]_] := Exp[-\[Eta]];
(*Implementation of Initial conditions*)
Do[f[[1, i + 1]] = f0[1, \[Eta]Grid[[i + 1]]];
F1[[1, i + 1]] = F10[1, \[Eta]Grid[[i + 1]]];
F2[[1, i + 1]] = F20[1, \[Eta]Grid[[i + 1]]];
G[[1, i + 1]] = G0[1, \[Eta]Grid[[i + 1]]];
H[[1, i + 1]] = H0[1, \[Eta]Grid[[i + 1]]];
S[[1, i + 1]] = S0[1, \[Eta]Grid[[i + 1]]];
P[[1, i + 1]] = P0[1, \[Eta]Grid[[i + 1]]];, {i, 0, N\[Eta]}];
(*Enforce boundary conditions*)
For[i = 1, i <= N\[Xi], i++,(*At \[Eta]=0*)
F1[[i, 1]] = 0.0; F2[[i, 1]] = 1.0; G[[i, 1]] = 1.0;
H[[i, 1]] = 1.0; S[[i, 1]] = 1.0; P[[i, 1]] = 1.0;
(*At \[Eta]=\[Eta]max*)
F1[[i, N\[Eta] + 1]] = 1.0; F2[[i, N\[Eta] + 1]] = 0.0;
G[[i, N\[Eta] + 1]] = 0.0; H[[i, N\[Eta] + 1]] = 0.0;
S[[i, N\[Eta] + 1]] = 0.0; P[[i, N\[Eta] + 1]] = 0.0;];
(*Now Define the Forward and Backward Derivatives Involved in \
Equations*)
f\[Xi][i_, j_] := If[i > 1, (f[[i, j]] - f[[i - 1, j]])/d\[Xi], 0];
F1\[Xi][i_, j_] :=
If[i > 1, (F1[[i, j]] - F1[[i - 1, j]])/d\[Xi], 0];
F2\[Xi][i_, j_] :=
If[i > 1, (F2[[i, j]] - F2[[i - 1, j]])/d\[Xi], 0];
G\[Xi][i_, j_] := If[i > 1, (G[[i, j]] - G[[i - 1, j]])/d\[Xi], 0];
H\[Xi][i_, j_] := If[i > 1, (H[[i, j]] - H[[i - 1, j]])/d\[Xi], 0];
S\[Xi][i_, j_] := If[i > 1, (S[[i, j]] - S[[i - 1, j]])/d\[Xi], 0];
P\[Xi][i_, j_] := If[i > 1, (P[[i, j]] - P[[i - 1, j]])/d\[Xi], 0];
P\[Eta][i_, j_] :=
If[1 < j < Length[\[Eta]Grid],
N[(P[[i, j + 1]] - P[[i, j - 1]])/(2*d\[Eta])], 0];
G\[Eta][i_, j_] :=
If[1 < j < Length[\[Eta]Grid],
N[(G[[i, j + 1]] - G[[i, j - 1]])/(2*d\[Eta])], 0];
G\[Eta]\[Eta][i_, j_] :=
If[1 < j < Length[\[Eta]Grid],
N[(G[[i, j + 1]] - 2*G[[i, j]] + G[[i, j - 1]])/(d\[Eta]^2)], 0];
f\[Xi]Grid =
Table[f\[Xi][i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
F1\[Xi]Grid =
Table[F1\[Xi][i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
F2\[Xi]Grid =
Table[F2\[Xi][i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
G\[Xi]Grid =
Table[G\[Xi][i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
H\[Xi]Grid =
Table[H\[Xi][i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
S\[Xi]Grid =
Table[S\[Xi][i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
P\[Xi]Grid =
Table[P\[Xi][i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
P\[Eta]Grid =
Table[P\[Eta][i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
G\[Eta]Grid =
Table[G\[Eta][i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
G\[Eta]\[Eta]Grid =
Table[G\[Eta]\[Eta][i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
Print["f\[Xi]Grid=", MatrixForm[f\[Xi]Grid[[All, 2 ;; -2]]]];
Print["F1\[Xi]Grid=", MatrixForm[F1\[Xi]Grid[[All, 2 ;; -2]]]];
Print["F2\[Xi]Grid=", MatrixForm[F2\[Xi]Grid[[All, 2 ;; -2]]]];
Print["G\[Xi]Grid=", MatrixForm[G\[Xi]Grid[[All, 2 ;; -2]]]];
Print["H\[Xi]Grid=", MatrixForm[H\[Xi]Grid[[All, 2 ;; -2]]]];
Print["S\[Xi]Grid=", MatrixForm[S\[Xi]Grid[[All, 2 ;; -2]]]];
Print["P\[Xi]Grid=", MatrixForm[P\[Xi]Grid[[All, 2 ;; -2]]]];
Print["P\[Eta]Grid=", MatrixForm[P\[Eta]Grid[[All, 2 ;; -2]]]];
Print["G\[Eta]Grid=", MatrixForm[G\[Eta]Grid[[All, 2 ;; -2]]]];
Print["G\[Eta]\[Eta]Grid=",
MatrixForm[G\[Eta]\[Eta]Grid[[All, 2 ;; -2]]]];
(*Now Define the Coefficients in Each Equation*)
A1[i_, j_] := f[[i, j]] + 2*B[[i]]*f\[Xi][i, j];
A2[i_, j_] := -(2 \[Beta][[i]]*F1[[i, j]] +
2 (B[[i]]*F1\[Xi][i, j]));
A3[i_, j_] := -2*B[[i]]*F1[[i, j]];
A4[i_, j_] := -2*\[Alpha][[i]]*F2[[i, j]];
A5[i_, j_] := \[Lambda]*s[[i]];
A6[i_, j_] := \[Lambda]*s[[i]]*Nc1;
A7[i_, j_] := \[Lambda]*s[[i]]*Nc2;
A8[i_, j_] := -\[Lambda]*s[[i]]*Nr;
U1[i_, j_] := -\[Beta][[i]]*(1 + F1[[i, j]]^2) -
2*B[[i]]*F1[[i, j]]*F1\[Xi][i, j] + \[Alpha][[i]]*F2[[i, j]]^2;
B1[i_, j_] := A1[i, j];
B2[i_, j_] := -2*\[Beta][[i]]*F1[[i, j]];
B3[i_, j_] := -2*B[[i]]*F1[[i, j]];
B4[i_, j_] := -2*\[Beta][[i]]*F2[[i, j]] - 2*B[[i]]*F2\[Xi][i, j];
U2[i_, j_] := -2*\[Beta][[i]]*F1[[i, j]]*F2[[i, j]] -
2*B[[i]]*F1[[i, j]]*F2\[Xi][i, j];
C1[i_, j_] := Pr*(A1[i, j] + Nb*P\[Eta][i, j] + 2*Nt*G\[Eta][i, j]);
C2[i_, j_] := -2*B[[i]]*Pr*F1[[i, j]];
C3[i_, j_] := -2*B[[i]]*Pr*G\[Xi][i, j];
C4[i_, j_] := Pr*Nb*G\[Eta][i, j];
U3[i_, j_] :=
Pr*Nt*G\[Eta][i, j]^2 - 2*B[[i]]*Pr*G\[Xi][i, j]*F1[[i, j]] +
Pr*Nb*G\[Eta][i, j]*P\[Eta][i, j];
D1[i_, j_] := Sc1*A1[i, j];
D2[i_, j_] := -2*B[[i]]*Sc1*F1[[i, j]];
D3[i_, j_] := -2*B[[i]]*Sc1*H\[Xi][i, j];
U4[i_, j_] := -2*B[[i]]*Sc1*H\[Xi][i, j]*F1[[i, j]];
E1[i_, j_] := Sc2*A1[i, j];
E2[i_, j_] := -2*B[[i]]*Sc2*F1[[i, j]];
E3[i_, j_] := -2*B[[i]]*Sc2*S\[Xi][i, j];
U5[i_, j_] := -2*B[[i]]*Sc2*S\[Xi][i, j]*F1[[i, j]];
I1[i_, j_] := Le*A1[i, j];
I2[i_, j_] := -2*B[[i]]*Le*F1[[i, j]];
I3[i_, j_] := -2*B[[i]]*Le*P\[Xi][i, j];
I4[i_, j_] := Nt/Nb;
U6[i_, j_] := -2*B[[i]]*Le*P\[Xi][i, j]*F1[[i, j]];
A1Grid =
Table[A1[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
A2Grid =
Table[A2[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
A3Grid =
Table[A3[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
A4Grid =
Table[A4[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
A5Grid =
Table[A5[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
A6Grid =
Table[A6[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
A7Grid =
Table[A7[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
A8Grid =
Table[A8[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
U1Grid =
Table[U1[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
B1Grid =
Table[B1[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
B2Grid =
Table[B2[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
B3Grid =
Table[B3[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
B4Grid =
Table[B4[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
U2Grid =
Table[U2[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
C1Grid =
Table[C1[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
C2Grid =
Table[C2[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
C3Grid =
Table[C3[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
C4Grid =
Table[C4[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
U3Grid =
Table[U3[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
D1Grid =
Table[D1[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
D2Grid =
Table[D2[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
D3Grid =
Table[D3[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
U4Grid =
Table[U4[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
E1Grid =
Table[E1[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
E2Grid =
Table[E2[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
E3Grid =
Table[E3[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
U5Grid =
Table[U5[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
I1Grid =
Table[I1[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
I2Grid =
Table[I2[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
I3Grid =
Table[I3[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
I4Grid =
Table[I4[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
U6Grid =
Table[U6[i, j], {i, 1, Length[\[Xi]Grid]}, {j, 1,
Length[\[Eta]Grid]}];
Print["A1Grid=", MatrixForm[A1Grid]]
Print["A2Grid=", MatrixForm[A2Grid]]
Print["A3Grid=", MatrixForm[A3Grid]]
Print["A4Grid=", MatrixForm[A4Grid]]
Print["A5Grid=", MatrixForm[A5Grid]]
Print["A6Grid=", MatrixForm[A6Grid]]
Print["A7Grid=", MatrixForm[A7Grid]]
Print["A8Grid=", MatrixForm[A8Grid]]
Print["U1Grid=", MatrixForm[U1Grid]]
Print["B1Grid=", MatrixForm[B1Grid]]
Print["B2Grid=", MatrixForm[B2Grid]]
Print["B3Grid=", MatrixForm[B3Grid]]
Print["B4Grid=", MatrixForm[B4Grid]]
Print["U2Grid=", MatrixForm[U2Grid]]
Print["C1Grid=", MatrixForm[C1Grid]]
Print["C2Grid=", MatrixForm[C2Grid]]
Print["C3Grid=", MatrixForm[C3Grid]]
Print["C4Grid=", MatrixForm[C4Grid]]
Print["U3Grid=", MatrixForm[U3Grid]]
Print["D1Grid=", MatrixForm[D1Grid]]
Print["D2Grid=", MatrixForm[D2Grid]]
Print["D3Grid=", MatrixForm[D3Grid]]
Print["U4Grid=", MatrixForm[U4Grid]]
Print["E1Grid=", MatrixForm[E1Grid]]
Print["E2Grid=", MatrixForm[E2Grid]]
Print["E3Grid=", MatrixForm[E3Grid]]
Print["U5Grid=", MatrixForm[U5Grid]]
Print["I1Grid=", MatrixForm[I1Grid]]
Print["I2Grid=", MatrixForm[I2Grid]]
Print["I3Grid=", MatrixForm[I3Grid]]
Print["I4Grid=", MatrixForm[I4Grid]]
Print["U6Grid=", MatrixForm[U6Grid]]
(*Tridiagonal matrix constructor with variable coefficients*)
TridiagonalMatrix[n_, lower_, diag_, upper_] :=
SparseArray[{Band[{1, 1}] -> diag, Band[{2, 1}] -> lower,
Band[{1, 2}] -> upper}, {n, n}];
errorPlots = {};
For[n = 1, n <= N\[Xi], n++,
F1Prev = F1[[n, 2 ;; -2]];
F2Prev = F2[[n, 2 ;; -2]];
GPrev = G[[n, 2 ;; -2]];
HPrev = H[[n, 2 ;; -2]];
SPrev = S[[n, 2 ;; -2]];
PPrev = P[[n, 2 ;; -2]];
Print["--- \[Xi] step: ", n, " ---"];
Print["F1Prev = ", F1Prev];
Print["F2Prev = ", F2Prev];
Print["GPrev = ", GPrev];
Print["HPrev = ", HPrev];
Print["SPrev = ", SPrev];
Print["PPrev = ", PPrev];
(*\[Xi]-derivatives (backward)*)
F1\[Xi]Grid[[n, 2 ;; -2]];
F2\[Xi]Grid[[n, 2 ;; -2]];
G\[Xi]Grid[[n, 2 ;; -2]];
H\[Xi]Grid[[n, 2 ;; -2]];
S\[Xi]Grid[[n, 2 ;; -2]];
P\[Xi]Grid[[n, 2 ;; -2]];
P\[Eta]Grid[[n, 2 ;; -2]];
G\[Eta]Grid[[n, 2 ;; -2]];
G\[Eta]\[Eta]Grid[[n, 2 ;; -2]];
Print["F1\[Xi]Grid=", MatrixForm[F1\[Xi]Grid[[n, 2 ;; -2]]]];
Print["F2\[Xi]Grid=", MatrixForm[F2\[Xi]Grid[[n, 2 ;; -2]]]];
Print["G\[Xi]Grid=", MatrixForm[G\[Xi]Grid[[n, 2 ;; -2]]]];
Print["H\[Xi]Grid=", MatrixForm[H\[Xi]Grid[[n, 2 ;; -2]]]];
Print["S\[Xi]Grid=", MatrixForm[S\[Xi]Grid[[n, 2 ;; -2]]]];
Print["P\[Xi]Grid=", MatrixForm[P\[Xi]Grid[[n, 2 ;; -2]]]];
Print["P\[Eta]Grid=", MatrixForm[P\[Eta]Grid[[n, 2 ;; -2]]]];
Print["G\[Eta]Grid=", MatrixForm[G\[Eta]Grid[[n, 2 ;; -2]]]];
Print["G\[Eta]\[Eta]Grid=",
MatrixForm[G\[Eta]\[Eta]Grid[[n, 2 ;; -2]]]];
(*Coefficients*)
A1coeff = Table[A1[n, j], {j, 2, N\[Eta]}];
A2coeff = Table[A2[n, j], {j, 2, N\[Eta]}];
A3coeff = Table[A3[n, j], {j, 2, N\[Eta]}];
A4coeff = Table[A4[n, j], {j, 2, N\[Eta]}];
A5coeff = Table[A5[n, j], {j, 2, N\[Eta]}];
A6coeff = Table[A6[n, j], {j, 2, N\[Eta]}];
A7coeff = Table[A7[n, j], {j, 2, N\[Eta]}];
A8coeff = Table[A8[n, j], {j, 2, N\[Eta]}];
U1coeff = Table[U1[n, j], {j, 2, N\[Eta]}];
B1coeff = Table[B1[n, j], {j, 2, N\[Eta]}];
B2coeff = Table[B2[n, j], {j, 2, N\[Eta]}];
B3coeff = Table[B3[n, j], {j, 2, N\[Eta]}];
B4coeff = Table[B4[n, j], {j, 2, N\[Eta]}];
U2coeff = Table[U2[n, j], {j, 2, N\[Eta]}];
C1coeff = Table[C1[n, j], {j, 2, N\[Eta]}];
C2coeff = Table[C2[n, j], {j, 2, N\[Eta]}];
C3coeff = Table[C3[n, j], {j, 2, N\[Eta]}];
C4coeff = Table[C4[n, j], {j, 2, N\[Eta]}];
U3coeff = Table[U3[n, j], {j, 2, N\[Eta]}];
D1coeff = Table[D1[n, j], {j, 2, N\[Eta]}];
D2coeff = Table[D2[n, j], {j, 2, N\[Eta]}];
D3coeff = Table[D3[n, j], {j, 2, N\[Eta]}];
U4coeff = Table[U4[n, j], {j, 2, N\[Eta]}];
E1coeff = Table[E1[n, j], {j, 2, N\[Eta]}];
E2coeff = Table[E2[n, j], {j, 2, N\[Eta]}];
E3coeff = Table[E3[n, j], {j, 2, N\[Eta]}];
U5coeff = Table[U5[n, j], {j, 2, N\[Eta]}];
I1coeff = Table[I1[n, j], {j, 2, N\[Eta]}];
I2coeff = Table[I2[n, j], {j, 2, N\[Eta]}];
I3coeff = Table[I3[n, j], {j, 2, N\[Eta]}];
I4coeff = Table[I4[n, j], {j, 2, N\[Eta]}];
U6coeff = Table[U6[n, j], {j, 2, N\[Eta]}];
(*System Matrices*)
AF1 = TridiagonalMatrix[N\[Eta] - 1,
Table[1 - A1coeff[[j]]*d\[Eta]/2, {j, 2, N\[Eta] - 2}],
Table[-2 + A2coeff[[j]]*d\[Eta]^2, {j, 2, N\[Eta] - 2}],
Table[1 + A1coeff[[j + 1]]*d\[Eta]/2, {j, 2, N\[Eta] - 2}]];
AF2 = TridiagonalMatrix[N\[Eta] - 1,
Table[1 - B1coeff[[j]]*d\[Eta]/2, {j, 2, N\[Eta] - 2}],
Table[-2 + B2coeff[[j]]*d\[Eta]^2, {j, 2, N\[Eta] - 2}],
Table[1 + B1coeff[[j + 1]]*d\[Eta]/2, {j, 2, N\[Eta] - 2}]];
AG = TridiagonalMatrix[N\[Eta] - 1,
Table[1 - C1coeff[[j]]*d\[Eta]/2, {j, 2, N\[Eta] - 2}],
Table[-2, {j, 2, N\[Eta] - 1}],
Table[1 + C1coeff[[j + 1]]*d\[Eta]/2, {j, 2, N\[Eta] - 2}]];
AH = TridiagonalMatrix[N\[Eta] - 1,
Table[1 - D1coeff[[j]]*d\[Eta]/2, {j, 2, N\[Eta] - 2}],
Table[-2, {j, 2, N\[Eta] - 1}],
Table[1 + D1coeff[[j + 1]]*d\[Eta]/2, {j, 2, N\[Eta] - 2}]];
AS = TridiagonalMatrix[N\[Eta] - 1,
Table[1 - E1coeff[[j]]*d\[Eta]/2, {j, 2, N\[Eta] - 2}],
Table[-2, {j, 2, N\[Eta] - 1}],
Table[1 + E1coeff[[j + 1]]*d\[Eta]/2, {j, 2, N\[Eta] - 2}]];
AP = TridiagonalMatrix[N\[Eta] - 1,
Table[1 - I1coeff[[j]]*d\[Eta]/2, {j, 2, N\[Eta] - 2}],
Table[-2, {j, 2, N\[Eta] - 1}],
Table[1 + I1coeff[[j + 1]]*d\[Eta]/2, {j, 2, N\[Eta] - 1}]];
(*Coupling Matrices*)
B12 = DiagonalMatrix[N\[Eta] - 1]*A4coeff*d\[Eta]^2;
B21 = DiagonalMatrix[N\[Eta] - 1]*B4coeff*d\[Eta]^2;
(*Right-Hand Sides*)
b1 = d\[Eta]^2*(U1coeff - A3coeff*F1\[Xi] - A5coeff*GPrev -
A6coeff*HPrev - A7coeff*SPrev - A8coeff*PPrev);
b2 = d\[Eta]^2*(U2coeff - B3coeff*F2\[Xi]);
b3 = d\[Eta]^2*(U3coeff - C2coeff*G\[Xi] - C3coeff*F1Prev -
C4coeff*G\[Eta]\[Eta]);
b4 = d\[Eta]^2*(U4coeff - D2coeff*H\[Xi] - D3coeff*F1Prev);
b5 = d\[Eta]^2*(U5coeff - E2coeff*S\[Xi] - E3coeff*F1Prev);
b6 = d\[Eta]^2*(U6coeff - I2coeff*P\[Xi] - I3coeff*F1Prev -
I4coeff*G\[Eta]\[Eta]);
(*Solve block system using Varga's method*){F1New, F2New, GNew, HNew,
SNew, PNew} =
VargaSolve[AF1, AF2, AG, AH, AS, AP, B12, B21, b1, b2, b3, b4, b5,
b6, F1Prev, F2Prev, GPrev, HPrev, SPrev, PPrev, 10^-4];
(*Update solution*)
F1[[n, 2 ;; -2]] = F1New;
F2[[n, 2 ;; -2]] = F2New;
G[[n, 2 ;; -2]] = GNew;
H[[n, 2 ;; -2]] = HNew;
S[[n, 2 ;; -2]] = SNew;
P[[n, 2 ;; -2]] = PNew;
F1[[n, All]] = Join[{0.0}, F1New, {1.0}];
F2[[n, All]] = Join[{1.0}, F2New, {0.0}];
G[[n, All]] = Join[{1.0}, GNew, {0.0}];
H[[n, All]] = Join[{1.0}, HNew, {0.0}];
S[[n, All]] = Join[{1.0}, SNew, {0.0}];
P[[n, All]] = Join[{1.0}, PNew, {0.0}];
(*Update f via trapezoidal rule again*)
For[j = 2, j <= Length[\[Eta]Grid], j++,
f[[n, j]] =
f[[n, 1]] +
Sum[0.5*(F1[[n, k]] + F1[[n, k - 1]])*d\[Eta], {k, 2, j}]];]
(*Error tracking*)
If[n == N\[Xi],
AppendTo[errorPlots,
ListLinePlot[errorList,
PlotLabel ->
"Error vs Iteration (Time step " <> ToString[n] <> ")",
Frame -> True, Axes -> True,
FrameLabel -> {"Iteration", "Max Residual Error"},
PlotStyle -> Red, ImageSize -> Medium]]];
(*Plot final solution*)
solPlot =
ListLinePlot[{F1[[-1]], F2[[-1]], G[[-1]], H[[-1]], S[[-1]],
P[[-1]]},
PlotLegends -> {"F1(\[Xi],\[Eta])", "F2(\[Xi],\[Eta])",
"G(\[Xi],\[Eta])", "H(\[Xi],\[Eta])", "S(\[Xi],\[Eta])",
"P(\[Xi],\[Eta])"}, PlotRange -> All,
PlotLabel -> "Final solution", Frame -> True, Axes -> True,
ImageSize -> Large];
(*Show both plots*)
Column[{solPlot, Sequence @@ errorPlots}]
I have a lot of Print command because I want to see where the problem is: But actually not so much expert of Mathematica
All the coefficients used in equations are defined in code structure.
