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

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?

• If you try a non-FEM method by expressing your boundary conditions as equations (i.e. bc = {L[0] == 1, L[100] == 0};) and removing the Method call inside NDSolve, it at least tries to evaluate. However, it ran for a few minutes on my computer without returning a result, probably because of the complexity of the problem, and I aborted it. Perhaps you can wait longer and see. Jun 18 '18 at 23:25
• @MarcoB Yeah, I'm running into the same problem where if I left it to run and nothing has happened after 1.5 hours. So, does this mean that the FEM method is not compatible with my problem? Is there a way to have it solve using another method? Jun 19 '18 at 4:37
• I don’t know enough about NDSolve to give you a better answer. There is, however, A LOT of information on this site on solving differential equations in general, and on FEM as well. I suggest a deep dive into some older answers first. Jun 19 '18 at 4:45
• A simple mistake is, your Plot won't work even if you successfully solve the equation, for more information, check document of NDSolve and ReplaceAll carefully. Then, NDSolve is very slow on this problem because the solution is severely oscillating, try e.g. sol = NDSolveValue[{eqn, L[0] == 1, L'[0] == 0.1}, L, {x, 0, 100}]; Plot[sol[x], {x, 0, 0.2}] Is it suppose to be like this? If not, you should double check the underlying model. Jun 19 '18 at 6:01
• This is a nonlinear stationary (i.e. not time dependent) PDE. This can not currently (Version 11.3) solved directly by NDSolve. Even if it were, one thing that is missing is an initial guess of L[x]. Also, are you sure the equations are correct? Jun 19 '18 at 8:51

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["NDSolveFEM"]

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 = NDSolveVariableData[{"DependentVariables", "Space"} -> {{L}, {x}}];
sd = NDSolveSolutionData[{"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"];

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];
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.

• Thank you! I modified my constants a bit and the r term. When I tried running it, it seems to freeze. Is there other parts in the code that I would need to modify if I modified any of the constants/equations? Jun 20 '18 at 20:40
• I found a potential cause for the freeze. Please try again. Jun 20 '18 at 20:43
• Unless I'm supposed to let it run for awhile, it still seems to be the same. Could the drastic change in the units for the constants and the way I re-expressed the r term may have affected it? Jun 20 '18 at 20:51
• Yes of course. Try to adapt them slowly. Jun 20 '18 at 20:54
• Okay. So it looks like changing my r term wasn't the problem. It was me changing my units from nM to M and micrometer^2 to meter^2. I'm a little confused following the code and knowing which part I should adapt based on the change in the constants. Jun 20 '18 at 21:02