# 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:
 Carl Woll 4 Votes 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}];