2
$\begingroup$

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

enter image description here

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?

$\endgroup$
9
  • $\begingroup$ 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. $\endgroup$
    – MarcoB
    Jun 18 '18 at 23:25
  • $\begingroup$ @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? $\endgroup$
    – BaiSango
    Jun 19 '18 at 4:37
  • $\begingroup$ 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. $\endgroup$
    – MarcoB
    Jun 19 '18 at 4:45
  • $\begingroup$ 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. $\endgroup$
    – xzczd
    Jun 19 '18 at 6:01
  • 2
    $\begingroup$ 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? $\endgroup$
    – user21
    Jun 19 '18 at 8:51
3
$\begingroup$

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}
 ]

enter image description here

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.

$\endgroup$
8
  • $\begingroup$ 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? $\endgroup$
    – BaiSango
    Jun 20 '18 at 20:40
  • $\begingroup$ I found a potential cause for the freeze. Please try again. $\endgroup$ Jun 20 '18 at 20:43
  • $\begingroup$ 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? $\endgroup$
    – BaiSango
    Jun 20 '18 at 20:51
  • $\begingroup$ Yes of course. Try to adapt them slowly. $\endgroup$ Jun 20 '18 at 20:54
  • $\begingroup$ 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. $\endgroup$
    – BaiSango
    Jun 20 '18 at 21:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.