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.
NDSolve
gives the error, the first or second? How isPlot
related to the problem? And so forth.) $\endgroup$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$