Kernel forever running with no error message. Solve NDSolve without FEM?

Question

I am trying to solve the following reaction-diffusion equation where R_L represents the reaction rate of L.

Clear["Global`*"]
k1 = 0.00193*10^-9;
k2 = 0.00255*10^-9;
k3 = 4.09*10^-9;
d1 = 0.00700;
d2 = 3.95*^-5;
d3 = 2.26;
K1 = 3.63*10^-9;
K2 = 0.155*10^-9;
K3 = 0.553*10^-9;
K4 = 9.01*10^-12;
Ki = 0.139;
NT = 1.70*10^-12;
L0 = 1*10^-9;
n = (200*^9)/(1.3*^-9);
NA = 6.02*^23;

r = -(K3*K4*
      Ki*(K2 + 
        L[x] - ((8*K2*NT*L[x]^2 + K3*K4*Ki*L[x]^2 + 8*K2*K3*NT*L[x] + 
             K2^2*K3*K4*Ki + 8*K2^2*K3*Ki*NT + 2*K2*K3*K4*Ki*L[x] + 
             8*K2*K3*Ki*NT*L[x])/(K3*K4*Ki))^(1/2)))/(4*(K3*L[x] + 
       L[x]^2 + K2*K3*Ki + K3*Ki*L[x]));

R = (L[x]/K1);
rr = (1/K4) r^2;
rR = (L[x]/((2*K2)*K4)) r^2;
pR = (L[x]/((2*K2)*K4*Ki)) r^2;
RR = (L[x]^2/(K2*K3*K4*Ki)) r^2;

r = (-(1 + (L[x]/
          K2)) + \[Sqrt]((1 + (L[x]/
             K2))^2 - (4*((2/
              K4)*((1 + (L[x]/
                  K2))*((1 + (1/
                    Ki))*((1 + (L[x]/K3))))))*(-NT))))/(2*((2/
         K4)*((1 + (L[x]/K2))*((1 + (1/Ki))*((1 + (L[x]/K3)))))));

myL = ((2*d1*(NT - r - rr - rR - pR - RR)) - (2*k1*L[x]*r)) + ((2*
      d2*(NT - r - R - rr - pR - RR)) - (k2*
      L[x]*(NT - r - R - rR - pR - RR))) + ((d3*(NT - r - R - rr - 
        rR - pR)) - (2*k3*L[x]*(NT - r - R - rr - rR - RR)))

(*---------------------------------------------------------------*)

diffCo = 1*10^-6;

bc = {L[0] == 1, L[100] == 0};

eqn = diffCo*L''[x] + myL == 0;

solNDSolve = 
 NDSolve[{eqn, bc}, L, {x, 0, 100}, 


  Plot[solNDSolve[x], {x, 0, 100}, PlotRange -> Full]

I originally used

 bc = {DirichletCondition[L[x] == 1, x == 0], 
      DirichletCondition[L[x] == 0, x == 100]};

which gave me back the errors:

CoefficientArrays::poly: -((0.00106107 L (-1-6.45161 L+Sqrt[12.3687 (1+Times[<<2>>]) (1+Times[<<2>>])+(1+Times[<<2>>])^2]))/((1+1.80832 L) (1+6.45161 L)))+2.26 (<<1>>)-0.00255 L (<<1>>)+0.014 (<<1>>)+0.000079 (<<1>>)-8.18 L (<<1>>)+L$45220/1000000 is not a polynomial.

and

NDSolve::femper: PDE parsing error of {-((0.00106107 L (-1-6.45161 L+Sqrt[12.3687 (1+Times[<<2>>]) (1+Times[<<2>>])+(1+Times[<<2>>])^2]))/((1+1.80832 L) (1+6.45161 L)))+<<7>>+L$45220/1000000}. Inconsistent equation dimensions.

when I switched it to

bc = {L[0] == 1, L[100] == 0};

I don't get the error message anymore, however, now it just runs without stopping/giving any error messages.

Solving using FEM doesn't seem to be working. Is there another way to solve this problem without using FEM?

Answer 1

In the meantime, I set up a simple Newton solver in order to solve these equations. It seems to work quite robustly as long as the initial values are between 0 and 1. Line search is crucial here in order to establish convergence.

Needs["NDSolve`FEM`"]

k1 = 0.00193;
k2 = 0.00255;
k3 = 4.09;
d1 = 0.00700;
d2 = 3.95*^-5;
d3 = 2.26;
K1 = 3.63;
K2 = 0.155;
K3 = 0.553;
K4 = 9.01;
Ki = 0.139;
NT = 1.70;
L0 = 1;
n = (200*^9)/(1.3*^-9);
NA = 6.02*^23;

R = (L[x]/K1);
rr = (1/K4) r^2;
rR = (L[x]/((2*K2)*K4)) r^2;
pR = (L[x]/((2*K2)*K4*Ki)) r^2;
RR = (L[x]^2/(K2*K3*K4*Ki)) r^2;

r = (-(1 + (L[x]/K2)) + √((1 + (L[x]/K2))^2 - (4*((2/K4)*((1 + (L[x]/K2))*((1 + (1/Ki))*((1 + (L[x]/K3))))))*(-NT))))/(2*((2/K4)*((1 + (L[x]/K2))*((1 + (1/Ki))*((1 + (L[x]/K3)))))));

myL = Simplify[((2*d1*(NT - r - rr - rR - pR - RR)) - (2*k1*L[x]* r)) + ((2*d2*(NT - r - R - rr - pR - RR)) - (k2* L[x]*(NT - r - R - rR - pR - RR))) + ((d3*(NT - r - R - rr - rR - pR)) - (2*k3*L[x]*(NT - r - R - rr - rR - RR)))];

reaction = L \[Function] Evaluate[Simplify[myL /. L[x] -> L]];
With[{diag = Simplify[D[reaction[L], L]]},
 Dreaction = Evaluate[L] \[Function] DiagonalMatrix[SparseArray[Flatten@diag]];
 ]

α = 1.;
β = 0.;
(*diffCo=1*10^-6;*)
diffCo = 10;
reg = ToElementMesh[DiscretizeRegion[Line[{{0}, {100}}]], "MeshOrder" -> 1, MaxCellMeasure -> .05];
bc = {
   DirichletCondition[L[x] == α, x == 0],
   DirichletCondition[L[x] == β, x == 100]
   };
vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{L}, {x}}];
sd = NDSolve`SolutionData[{"Space"} -> {reg}];
cdata = InitializePDECoefficients[vd, sd,
   "DiffusionCoefficients" -> {{diffCo}},
   "MassCoefficients" -> {{1}}
];
bcdata = InitializeBoundaryConditions[vd, sd, bc];
mdata = InitializePDEMethodData[vd, sd];

(*Discretization*)
dpde = DiscretizePDE[cdata, mdata, sd];
dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd];
{load, stiffness, damping, mass} = dpde["All"];
DeployBoundaryConditions[{load, stiffness}, dbc];

F = u \[Function] Module[{b},
    b = Partition[stiffness.u + mass.reaction[u], 1];
    DeployBoundaryConditions[{b, stiffness}, dbc];
    b[[1]] = u[[1]] - α;
    b[[-1]] = u[[-1]] - β;
    Flatten[b]
    ];

F' = u \[Function] Module[{A}, A = stiffness + mass.Dreaction[u];
    DeployBoundaryConditions[{load, A}, dbc];
    A
    ];

NewtonStep = u \[Function] Module[{F0, Ft, θ0, θt, δu, t, γ, ut, bool, σ},
    γ = 0.1;
    σ = 0.01;

    F0 = F[u];
    θ0 = F0.(mass.F0);
    δu = -LinearSolve[F'[u], F0, Method -> "Banded"];

    t = 1.;
    ut = u + t δu;
    Ft = F[ut];
    θt = Ft.(mass.Ft);
    bool = Not[θt ∈ Reals];
    If[! bool, bool = θt >= (1. - σ t)];
    (* Line search*)
    While[bool,
     t *= γ;
     ut = u + t δu;
     Ft = F[ut];
     θt = Ft.(mass.Ft);
     bool = Not[θt ∈ Reals];
     If[! bool, bool = θt >= (1. - σ t)];
     ];
    ut
    ];

Setting up an initial condition and let Newton's method run.

u = RandomReal[{0., 1.}, {Length[mass]}];
u[[1]] = α;
u[[-1]] = β;
data = FixedPointList[NewtonStep, u, SameTest -> (Max[Abs[#1 - #2]] < 1*^-14 &)];

Visualization of the solving process:

Manipulate[
 ListLinePlot[Flatten[data[[i]]], PlotRange -> {0, 1}, 
  AxesLabel -> {"x", "L"}, PlotLabel -> "Step " <> ToString[i - 1]],
 {i, 1, Length[data], 1}
 ]

Note that I cranked up the diffusivity considerably. Otherwise the flanks are really steep (and for very small diffusivities, a lot of oscillaction is happening that should not be there) and one would hardly see anything.