(Subscript[SP, L]) equations

Interior equations

Expansion of solution in zonal spherical harmonics (1D Legendre polynomials)

ψ[μ_,L_]:=(∑_(n=0)^L)((2n+1)/2)LegendreP[n,μ] (Subscript[ϕ, n])[z]

L=3;

(Subscript[P, n])equations (replace Q by d/dz)

P[n_]=((n+1)/(2n+1))Q (Subscript[ϕ, n+1])[z] +(n/(2n+1))Q (Subscript[ϕ, n-1])[z]+(Subscript[Σ, a,n])[z](Subscript[ϕ, n])[z]==(Subscript[S, n])[z]

eqns=Table[P[n]/.If[n≥1, (Subscript[S, n])[z]->0,(Subscript[S, n])[z]->(Subscript[S, n])[z]],{n,0,L}]/.(Subscript[ϕ, L+1])[z]->0;

MatrixForm[eqns]

Expansion moments (odd/even)

oddn=Table[2n-1,{n,1,(L+1)/2}];
evenn=Table[2n,{n,1,(L+1)/2}];

oddmom=Table[(Subscript[ϕ, n])[z],{n,0,L}][[evenn]]
evenmom=Table[(Subscript[ϕ, n])[z],{n,0,L}][[oddn]]

Reduction to the “simplified (Subscript[P, L]) form” (second order eqns.)

express odd moments using the even ones

oddWrtEven=Solve[eqns[[evenn]],oddmom]

insert into the equations

SPL=eqns[[oddn]]/.oddWrtEven[[1]]/.(Subscript[S, 0])[z]->(Subscript[νΣ, f])(Subscript[ϕ, 0])[z];

MatrixForm[SPL]

Definition of “pseudo-fluxes” (Subscript[F, i]) and “pseudo-diffusion-coefficients” (Subscript[D, i])
such that “pseudo-currents” may be defined according to “pseudo-Fick’s-law” as (Subscript[D, 2i+1]) d(Subscript[F, i])/dz (here (Subscript[D, i])Q (Subscript[F, i]))

evenPseudoFluxes=Table[(n+1) (Subscript[ϕ, n])[z] +(n+2) (Subscript[ϕ, n+2])[z]==(Subscript[F, n])[z],{n,0,L,2}]/.(Subscript[ϕ, L+1])[z]->0
pseudoDiffusions=Table[(1/((2n+1)(Subscript[Σ, a,n])[z]))->(Subscript[D, n])[z],{n,1,L,2}]

Even expansion moments expressed using the even pseudo-fluxes
the zero-th moment is the scalar flux

Table[(Subscript[ϕ, 2i])[z]/.Solve[Eliminate[evenPseudoFluxes,Drop[evenmom,{i+1}]],(Subscript[ϕ, 2i])[z]]//First//Expand,{i,0,((L-1)/2)}]//MatrixForm

Expressing the (Subscript[SP, N])equations in the operator form A {(Subscript[F, i])}={(Subscript[b, i])}, where the streaming matrix A=diag((d/(dz))(Subscript[D, 2i+1])(d/(dz)))

lower=Table[(n/(2n+1)),{n,2,L,2}]
centr=Table[((n+1)/(2n+1)),{n,0,L,2}]

Adata=If[L>1,{Band[{2,1}]->lower,Band[{1,1}]->centr},{Band[{1,1}]->centr}];

A=SparseArray[Adata,Floor[L/2]+1];

A//Normal//MatrixForm

rhs=Table[(Subscript[Σ, a,n])[z](Subscript[ϕ, n])[z],{n,0,L,2}];

Eliminating the even flux  expansion moments in favor of the unknowns: the even pseudo-fluxes (Subscript[F, i])

DelDDelFs=Table[(Subscript[DelDDelF, n]),{n,0,L,2}];

eqnsForDelDDelFs1=Thread[DelDDelFs==LinearSolve[A,rhs]]

eqnsForDelDDelFs2=Map[Join[{#},evenPseudoFluxes]&,eqnsForDelDDelFs1]

eqnsForDelDDelFs3=Map[Eliminate[#,evenmom]&,eqnsForDelDDelFs2]

coefAtDs=Thread[Coefficient[eqnsForDelDDelFs3/.Equal->Plus,DelDDelFs]]

eqnsForDelDDelFs4=Thread[Thread[Expand[eqnsForDelDDelFs3[[All,1]]/coefAtDs]]==Thread[Expand[eqnsForDelDDelFs3[[All,2]]/coefAtDs]]]

eqnsForDelDDelFs5=eqnsForDelDDelFs4/.Thread[DelDDelFs->Table[(Subscript[DelDDel, n])[z](Subscript[F, n])[z],{n,0,L,2}]]

Rewriting the equations into a form  m {(Subscript[F, i])}-{(Subscript[b, i])}==0,  where now the matrix m contains both the streaming and reaction operators

{b,m}=CoefficientArrays[eqnsForDelDDelFs5,Table[(Subscript[F, n])[z],{n,0,L,2}]]

Normal[b]//MatrixForm

Normal[m]//MatrixForm

Some formal changes to improve readability

F[A_]:=TraditionalForm[A/.{(Subscript[DelDDel, n_])[z]->(Subscript[DelDDel, n]),(Subscript[Σ, a,n__])[z]->(Subscript[Σ, n])}]

F[Normal[m]//MatrixForm]

Make sure that in the monoenergetic case, m is symmetric

Normal[m-Transpose[m]]

m2=-Normal[m]/.{(Subscript[DelDDel, n_])[z]->(Subscript[DelDDel, n]),(Subscript[Σ, a,n__])[z]->(Subscript[Σ, n])}

m2//MatrixForm

MatrixForm[m2/.(Subscript[DelDDel, _])->0]

Sm=0;

For[i=0,i<L,i+=2,
{b,s0}=CoefficientArrays[m2, (Subscript[Σ, i])];
c[i]=Flatten[ Normal[s0],{1,3}];
Print[MatrixForm[c[i]]];
Print[Eigenvalues[c[i]]];
Sm=Sm+c[i]
]

Eigenvalues[Sm]

MatrixForm[Sm]

Boundary conditions

Legendre expansion of the outgoing partial current

(SuperPlus[j])=Table[(∫_0^1)LegendreP[n,μ]ψ[μ,L]\[DifferentialD]μ,{n,1,L,2}]

expressed using only the even order moments (as in the simplified form of the (Subscript[P, n]) equations)

(SuperPlus[j])=(SuperPlus[j])/.oddWrtEven[[1]]

Extract the equations for pseudo-currents

DDelFs=Table[(Subscript[DDelF, n]),{n,0,L,2}]

eqnsForDDelFs1=Thread[DDelFs==2Map[Delete[#,-1]&,(SuperPlus[j])]]

Express the boundary equations in terms of only the even expansion moments and pseudo-currents

eqnsForDDelFs2=Map[Join[{#},evenPseudoFluxes]&,eqnsForDDelFs1]

The final Robin-type boundary conditions

bc1p=Map[Eliminate[#,evenmom]&,eqnsForDDelFs2]

bc2p=Solve[bc1p,DDelFs]//Flatten//Expand

The same as above, only for the incoming partial currents

(SuperMinus[j])=Table[(∫_(-1)^0)LegendreP[n,Abs[μ]]ψ[μ,L]\[DifferentialD]μ//Distribute,{n,1,L,2}]

(SuperMinus[j])=(SuperMinus[j])/.oddWrtEven[[1]]

eqnsForDDelFs1=Thread[DDelFs==-2Map[Delete[#,-1]&,(SuperMinus[j])]]

eqnsForDDelFs2=Map[Join[{#},evenPseudoFluxes]&,eqnsForDDelFs1]

bc1m=Map[Eliminate[#,evenmom]&,eqnsForDDelFs2]

bc2m=Solve[bc1m,DDelFs]//Flatten//Expand

curs = -DDelFs/.bc2m

ca=CoefficientArrays[curs,evenmom/.ϕ->F]

Normal[ca[[2]]]

N[Eigenvalues[ca[[2]]]]

Albedo boundary conditions

alpha=Solve[γ== ((1-α)/(2(1+α))),α][[1]]

aj=Thread[(SuperMinus[j])==α (SuperPlus[j])]

DDelFsPos=1+((L+1)/2);

DDelFsRep=Flatten[Table[{{{i,1,DDelFsPos}}->((DDelFs[[i]])/2),{{i,2,2,DDelFsPos}}->-((DDelFs[[i]])/2)},{i,1,((L+1)/2)}]]

aj2=ReplacePart[aj,DDelFsRep]

eqnsForDDelFs3=Map[Join[{#},evenPseudoFluxes]&,aj2]

bc1a=Map[Eliminate[#,evenmom]&,eqnsForDDelFs3]

bc2a=Solve[bc1a,DDelFs]//Flatten

bc3a=bc2a/.alpha//FullSimplify

cursa = -DDelFs/.bc3a

ca=CoefficientArrays[cursa,evenmom/.ϕ->F]

Normal[ca[[2]]]//MatrixForm

Normal[ca[[2]]/γ]//MatrixForm

Eigenvalues[Normal[ca[[2]]]]//N

γf[α_]:=((1-α)/(2(1+α)));
Plot[γf[α],{α,0,1}]