Message Boards Message Boards

Correct a NDSolve approach when an argument contains an InverseFunction?

GROUPS:

We have a thin crystal of thickness d illuminated uniformly from the left at an intensity lgti[0, t] of unity. The crystal is composed of a photo reactive species a which absorbs light. Inside the crystal at location x and time t the photo chemical reaction leads to a local concentration of a given as a[x, t] in the code below.

Code to generate the concentration, a[x, t], and light intensity, lgti[x, t], within the crystal is straightforward and is shown.

The cause of the errors generated using the code is that in some cases a[x, t] does not return a number.These exceptional cases are near a known limit of the InverseFunction and a[x, t] in these cases could be given the values unity if found. My problem is that using NumberQ inside a Module definition of a[x, t] always gives false because ligti[x, t] is unevaluated.

The starting code is shown but see the attached notebook.

ClearAll[a, x, t, eqns, \[Sigma]N, lgti, soln, t0, d, tmax]

\[Sigma]N = 1.78*10^4; d = 0.001; tmax = 3000.0;

a[x_, t_] := InverseFunction[(1.8996253051848473`/lgti[x, t] *
      (162.99559471365637` Log[1 + 8.98898678414097` (1 - #1)] - 
        172.98458149779734` Log[#1]) ) &
                                                            ] [t]

eqns = {D[lgti[x, t], x ] == - \[Sigma]N a[x, t ] lgti[x, t],(*    Beer's Law *)
                   lgti[0, t] == 1, 
   lgti[x, 0] == Exp[-\[Sigma]N x]
                    };

t0 = AbsoluteTime[];

soln = NDSolve[eqns, lgti, {x, 0, d}, {t, 0, tmax}, 
              MaxStepFraction -> 0.01]   [[1]];

Print[ToString[(AbsoluteTime[] - t0)/60] <> "   minutes"]

Any advice on how to code a[x, t] so that lgti[x, t] appears as a number within the body of the code would be welcome.

An alternate approach would also be well received.

Attachments:
POSTED BY: Mervin Hanson
Answer
1 year ago

I wouldn't use InverseFunction in the body of the PDE in this way, since NDSolve is perfectly capable of computing the inverse while solving the PDE. To make this post self-contained, I include a few of the OP details:

\[Sigma]N = 1.78*10^4; d = 0.001; tmax = 3000.0;

inversePureFunction = (1.8996253051848473`/lgti[x,t]*(162.99559471365637` Log[1+8.98898678414097` (1-#1)]-172.98458149779734` Log[#1]))&;

pde1 = D[lgti[x, t], x ] == - \[Sigma]N a[x, t ] lgti[x, t];

initialConditions = {lgti[0,t] == 1, lgti[x, 0] == Exp[- \[Sigma]N x]};

Now, the first step is to convert the InverseFunction into an implicit equation. Here is the implicit equation equivalent:

In[50]:= implicitEquation = Thread[
    lgti[x,t] (inversePureFunction[a[x,t]]==t),
    Equal
]
Out[50]= 1.89963 (162.996 Log[1+8.98899 (1-a[x,t])]-172.985 Log[a[x,t]])==t lgti[x,t]

The idea is to add a pde for a[x,t] to the NDSolve call. So, we need to come up with a pde for a[x,t] and initial conditions. Using pde1 and implicitEquation, we can solve for D[a[x,t],x]:

In[51]:= asol=Quiet[
    Solve[{D[implicitEquation,x],pde1},D[a[x,t],x],D[lgti[x,t],x]],
    Solve::bdomv
];
asol//InputForm
Out[52]//InputForm=
{{Derivative[1, 0][a][x, t] -> (938.0600000000017*t*a[x, t]^2*(-1.111247243322715 + 1.*a[x, t])*lgti[x, t])/(-19.244057828963523 + 1.*a[x, t])}}

where I use InputForm to avoid ambiguity (and Quiet the Solve::bdomv message). This results in the following pde for a[x,t]:

In[53]:= pde2=Equal@@asol[[1,1]];
pde2//InputForm
Out[54]//InputForm=
Derivative[1, 0][a][x, t] == (938.0600000000017*t*a[x, t]^2*(-1.111247243322715 + 1.*a[x, t])*lgti[x, t])/(-19.244057828963523 + 1.*a[x, t])

Now, we need initial conditions. When t==0, we have:

In[55]:= implicitEquation/.t->0
Out[55]= 1.89963 (162.996 Log[1+8.98899 (1-a[x,0])]-172.985 Log[a[x,0]])==0

which obviously has the solution:

a[x,0]==1

(although getting WL to prove this isn't easy). When x==0, we can just use the original InverseFunction, since it has no problems when x is 0:

In[56]:= a0t = InverseFunction[inversePureFunction/.x->0/.lgti[0,t]->1][t]
Out[56]= InverseFunction[1.89963 (162.996 Log[1+8.98899 (1-#1)]-172.985 Log[#1])&][t]

Finally, we can use NDSolve:

NDSolve[{pde1,pde2,lgti[x,0]==Exp[-\[Sigma]N x],lgti[0,t]==1,a[x,0]==1,a[0,t]==a0t},{a,lgti},{x,0,d},{t,0,tmax}];
POSTED BY: Carl Woll
Answer
11 months ago

Thank you

POSTED BY: Mervin Hanson
Answer
8 months ago

Thank you. I had given up hope. Unfortunately I will be tied up with a medical condition for a couple of months. With my current meds it takes me15-20 sec to find the time 09:37 + 3:30

Something to look forward to.

POSTED BY: Mervin Hanson
Answer
11 months ago

Group Abstract Group Abstract