0
$\begingroup$

I'm trying to solve a set of differential-algebraic equations with boundary conditions:

$A\frac{d^2\psi}{dx^2}=[c_{F}-2c_{D}(x)+c_{C}(x)]F$

$c_{D}(x)=0.1\exp(-0.2253 - 76.5596\psi(x) - 1.1149)$

$c_{C}(x)=0.2\exp(-1.3305 + 38.2798\psi(x) - 0.3459)$

$\psi\,'(0)=0,\psi(\frac{1}{2})=0$

With the following code, shooting method is used to find the Dirichlet BC at $x=0$ to use NDsolve to solve the DAE

Clear["Global`*"]

A = 5.31251*10^-10;
F = 96485.33252;
cF = 1000;

(*Obtain Dirichlet condition of psi at x=0 by shooting method*)
fun[y_?NumericQ] := Module[{eqns, sol, ics, x, psi, cD, cC},
   eqns = {A*D[psi[x], x, x] == (cF - 2*cD[x] + cC[x])*F,
     cD[x] == 0.1*Exp[-0.2253 - 76.5596*psi[x] - 1.1149],
     cC[x] == 0.2*Exp[-1.3305 + 38.2798*psi[x] - 0.3459]};
   ics = {psi[0] == y, psi'[x] == 0 /. x -> 0};
   sol = 
    NDSolve[{eqns, ics}, {psi[x], cD[x], cC[x]}, {x, 0, 1/2}][[1]];
   psi[x] /. sol /. x -> 1/2];
ic = FindRoot[fun[x], {x, 0}]

dae = {A*D[psi[x], x, x] == (cF - 2*cD[x] + cC[x])*F,
   cD[x] == 0.1*Exp[-0.2253 - 76.5596*psi[x] - 1.1149],
   cC[x] == 0.2*Exp[-1.3305 + 38.2798*psi[x] - 0.3459]};
bcs = {psi[0] == x /. ic[[1]], D[psi[x], x] == 0 /. x -> 0};

(*Solve the differential equation using NDSolve*)
sol = NDSolve[{dae, bcs}, {psi, cD, cC}, {x, 0, 1/2}];

The NDSolve returns the following error:

Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.

Even if the system of DAE is converted to an ODE to eliminate the need for shooting method, applying NDSolve gives the error:

Clear["Global`*"]

A = 5.31251*10^-10;
F = 96485.33252;
cF = 1000;

ode = {A*
     D[psi[x], x, 
      x] == (cF - 2*0.1*Exp[-0.2253 - 76.5596*psi[x] - 1.1149] + 
       0.2*Exp[-1.3305 + 38.2798*psi[x] - 0.3459])*F};
bcs = {psi[1/2] == 0, D[psi[x], x] == 0 /. x -> 0};

(*Solve the differential equation using NDSolve*)
sol = NDSolve[{ode, bcs}, {psi}, {x, 0, 1/2}];

NDSolve::ndsz: At x == 2.012990504102674`*^-9, step size is effectively zero; singularity or stiff system suspected.

Are there any methods that I can use to get around this? I'm new to numerical method, any help is greatly appreciated, thanks a lot.

$\endgroup$
4
  • $\begingroup$ You could try giving initial conditions for both values and derivatives of the functions. It requires an understanding of the problem, which I don't have at present. Knowing what gave rise to the system your solving might help. $\endgroup$
    – Michael E2
    Commented Feb 21, 2023 at 14:00
  • $\begingroup$ It also would be nice to isolate the problem for site users. (Which NDSolve gives the error, the first or second? How is Plot related to the problem? And so forth.) $\endgroup$
    – Michael E2
    Commented Feb 21, 2023 at 14:08
  • $\begingroup$ @MichaelE2 The first NDSolve failed, and even if I made the DAEs into an ODE (to eliminate the need for shooting method), NDSolve also failed to solve using the two BCs stated above, giving the error NDSolve::ndsz: At x == 2.012990504102674*^-9, step size is effectively zero; singularity or stiff system suspected. Plot is just for plotting the solution, in this case it isn't related to the problem of NDSolve. Is there any way to get around this? $\endgroup$
    – Jack L
    Commented Feb 22, 2023 at 5:05
  • $\begingroup$ @MichaelE2 This system arose from the calculation of concentration of ions in a membrane. Currently I only have the 2 BCs stated to work with, are there any methods I can implement for NDSolve for this? $\endgroup$
    – Jack L
    Commented Feb 22, 2023 at 5:14

1 Answer 1

0
$\begingroup$

Eliminating the unnecessary variables converts the first system from a DAE to an ODE, which gets rid of the problem complained of in the OP:

fun[y_?NumericQ] := 
  Module[{eqns, sol, ics, x, psi, cD, cC}, 
   eqns = {A*
       D[psi[x], x, 
        x] == (cF - 
         2*(0.1*Exp[-0.2253 - 76.5596*psi[x] - 1.1149]) + (0.2*
           Exp[-1.3305 + 38.2798*psi[x] - 0.3459]))*F};
   ics = {psi[0] == y, psi'[x] == 0 /. x -> 0};
   NDSolveValue[{eqns, ics}, {psi[x]}, {x, 0, 1/2}]
   ];

However, the second derivative of psi grows so fast, other problems arise, naturally, such as singularities.

The second derivative as a function of psi:

Plot[(cF - 
    2*(0.1*Exp[-0.2253 - 76.5596*psi[x] - 1.1149]) + (0.2*
      Exp[-1.3305 + 38.2798*psi[x] - 0.3459]))*F/A, {psi[x], 0, 1}]
$\endgroup$
3
  • $\begingroup$ With the function fun presented, it will be called by FindRoot in the subsequent line, and from there NDSolveVelue will fail, giving the error NDSolveValue::ndsz: At x$7057 == 2.012990266193402*^-9, step size is effectively zero; singularity or stiff system suspected. As mentioned by you, singularities may occur, hence do you mind suggesting any inbuilt numerical methods that can solve this problem? $\endgroup$
    – Jack L
    Commented Feb 22, 2023 at 5:15
  • $\begingroup$ @JackL The two plots show a large second derivative in the neighborhood of psi == 0 (right), which is a BC in your setup. The left plot shows rapid growth of the second derivative, suggesting that a singularity might be inevitable. It's a feature of the equations that any accurate numerical method would show. If there is a way to rescale psi, it might turn out that the singularity is a numerical artifact of poor scaling. Perhaps recheck the equations. $\endgroup$
    – Michael E2
    Commented Feb 22, 2023 at 11:54
  • $\begingroup$ Alright, thanks a lot for the suggestion, I'll look into it. $\endgroup$
    – Jack L
    Commented Feb 22, 2023 at 14:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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