# NDSolve::bcsol?

I'm writing a code to numerically solve a differential equation. It has to repeat the process several times since there is a variable involved which is unknown, and the equation is fit to some experimental data (just know that the equation has to solve many times). However, I'm getting an error:

NDSolve::bcsol: Could not solve for equations at boundary points from the boundary conditions.

The Mathematica help section says:

"This message is generated when the partial differential equation solver is unable to solve algebraically for the boundary conditions needed to compute a solution."

The code itself requires importing some data from excel, but here are the basics:

Initial data:

nldiss = 1.91;(*concentration of dissolved species in the pore fluid,
in mol/L*)
rhow0 = 1000;(*Density of water kg/m^2*)
R = 8.314; (*J/(K mol)*)
Temp = 293;(*Kelvin*)
Mw = 18/1000;(*g/mol*)
μw0 = -237000; L = 2.01/2(*mm*);
phi = 0.445;
diinitial = .05(*kg/(mm hr*);
wcm =0.5;
Vpaste = 1.7*10^-5;
nlw = 55;(*molar concentration of pure water*)
RHmax = (nlw/(nldiss + nlw));(*max RH in at S=1*)
alpha=0.548994;
beta=0.197269;
dilist050={0.37};
MaxTime=6;
tdata={0.086113, 0.175003, 0.263894, 0.352784, 0.441675, 0.530566, \
0.619456, 0.708347, 0.797238, 0.886128, 0.975019, 1.06391, 1.1528, \
1.2417, 1.33061, 1.41951, 1.5084, 1.5973, 1.68619, 1.77509, 1.86398, \
1.95288, 2.04177, 2.13066, 2.21956, 2.30845, 2.39735, 2.48624, \
2.57513, 2.66403, 2.75292, 2.84182, 2.93071, 3.0196, 3.1085, 3.19739, \
3.28629, 3.37518, 3.46407, 3.55297, 3.64186, 3.73076, 3.81965, \
3.90854, 3.99744, 4.08633, 4.17523, 4.26412, 4.35301, 4.44191, \
4.5308, 4.6197, 4.70859, 4.79748, 4.88638, 4.97527, 5.06417, 5.15306, \
5.24195, 5.33085, 5.41974, 5.50864, 5.59753, 5.68642, 5.77532, \
5.86421, 5.95311, 6.042, 6.13089, 6.21979, 6.30868, 6.39758, 6.48647, \
6.57536, 6.66426};

(*Boundary Conditons:*)

RHbound = RHmax - (RHmax - .85)*Tanh[t*100];
Sbound[t_] = Exp[-alpha*Log[(1 - Log[RHbound/RHmax]/beta)]];
Sbounddot[t_] = D[Sbound[t], t];
μbound[t_] = μw0 +
R*Temp*(Log[RHbound/1] + Log[1 + nldiss*(1/nlw)/Sbound[t]])/Mw;
interfacesboundmax =
NIntegrate[(μbound[tp] - μw0) Sbounddot[tp]/Sbound[tp], {tp,
0, MaxTime}];
ff[t_] = interfacesboundmax*Tanh[t*100];
IC1 = μ[0, x] == -237000;
BC1 = μ[t, 0] == μw0 +
R*Temp*(Log[RHbound/1] + Log[1 + nldiss*(1/nlw)/Sbound[t]])/Mw -
ff[t];
BC2 = μ[t, 0] == μw0 +
R*Temp*(Log[RHbound/1] + Log[1 + nldiss*(1/nlw)/Sbound[t]])/Mw -
ff[t];

Clear[μ, RH, AAA]

(*Set up differential equation:*)

Sμ = 1 - 0.011672402233208307*(-232380 - μ[t, x])^(1/3);
Sμdot = Simplify[D[Sμ, t]];
AAA = -(μ[t, x] - μw0) Sμdot/Sμ;
eq1[di_] = D[μ[t, x], t] == AAA + di Derivative[0, 2][μ][t, x];

(*Solve differential equation:*)

Clear[di, RHfit, RHavg, Savg, Saturation, Massloss, modeldata,
residuals, residualssquare, Rsquare]

fun[data_, deq_] := Module[{sol2, residuals},
sol2 = NDSolve[{deq, IC1, BC1, BC2}, μ, {t, 0, MaxTime}, {x, 0, 2 L}];
Saturation[t_, x_] =
1 - 0.011672402233208307*(-232380 -
Evaluate[μ[t, x] /. sol2])^(1/3);
Savg[t_] := 1/L Integrate[Saturation[t, x], {x, 0, L}];
Massloss[t_] := (Savg[t] - Savg) rhow0 phi Vpaste*1000;
modeldata = Map[Evaluate[Massloss], tdata][[All, 1]];
residuals = N[data[[All, 2]] - modeldata];
residualssquare = Total[residuals^2]
];

fun1[di_?NumberQ] := fun[datalist, eq1[di]]

Clear[dii, res]

{res, dii} =
Reap[FindMinimum[{fun1[di], .001 < di < 1}, {di, Last[dilist050]},
MaxIterations -> 8, AccuracyGoal -> 4, PrecisionGoal -> 3,
WorkingPrecision -> MachinePrecision, StepMonitor :> Sow[di],
Method -> InteriorPoint]]


Ok, so all of that code requires many variable inputs (such as R and Temp), but hopefully something in there might give a clue.

Here's the big issue: it spits out the NDSolve::bscol error, but this is not the first time that this code has been run, and it runs the exact same code earlier in my program with no problems (I have to run this piece of code 5 times in my program). I've checked the number values for boundary conditions and initial conditions to make sure they agree, and they do. I have no idea what this error is telling me except that it's a boundary condition problem. Does anyone know more information about this specific error?

• I tried to fix some of the formatting. Please check I didn't mistakenly change the code. -- The trouble at this point is that there are so many undefined terms, it is hard to know why the boundary conditions cannot be solved. (1) Even if it's "gross", copy/paste/execute in Mathematica is not hard to do. Or (2) perhaps you could post what IC1, BC1, BC2 each evaluates to. – Michael E2 May 15 '15 at 16:40
• It appears that BC1, BC2 are the same. Are they meant to be? In any case, I think that's the source of the error message. – Michael E2 May 15 '15 at 17:53
• This yields a solution without errors: NDSolve[{eq1[Last[dilist050]], Derivative[0, 1][\[Mu]][t, 0] == 0, IC1, BC1}, \[Mu], {t, 0, MaxTime}, {x, 0, 2*L}]. Of course the BC Derivative[0, 1][\[Mu]][t, 0] == 0 was just made up. (I don't think there's any algebraic difference between X == 0 and 0.999 X == 0, if that's what you meant.) – Michael E2 May 15 '15 at 18:04