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[0]) 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?
Thanks in advance! Josh
IC1, BC1, BC2
each evaluates to. $\endgroup$BC1, BC2
are the same. Are they meant to be? In any case, I think that's the source of the error message. $\endgroup$NDSolve[{eq1[Last[dilist050]], Derivative[0, 1][\[Mu]][t, 0] == 0, IC1, BC1}, \[Mu], {t, 0, MaxTime}, {x, 0, 2*L}]
. Of course the BCDerivative[0, 1][\[Mu]][t, 0] == 0
was just made up. (I don't think there's any algebraic difference betweenX == 0
and0.999 X == 0
, if that's what you meant.) $\endgroup$