0
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Michael E2 May 15 '15 at 16:40
  • $\begingroup$ Ok, I think I have added enough data to make it "run", i.e. it still doesn't work, but it's giving the error that I get. $\endgroup$ – Hogy88 May 15 '15 at 17:37
  • $\begingroup$ 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. $\endgroup$ – Michael E2 May 15 '15 at 17:53
  • $\begingroup$ Yes, BC1 and BC2 are the same. I just tried changing BC1 to 0.999*BC1 (so just a minor change) and it gave the same warning (plus a warning stating that the boundary condition and initial condition are inconsistent ) $\endgroup$ – Hogy88 May 15 '15 at 17:58
  • 1
    $\begingroup$ 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.) $\endgroup$ – Michael E2 May 15 '15 at 18:04
2
$\begingroup$

The solution to this problem was a simple one: BC1 was defined at [t,0], and BC2 was also defined at [t,0]. Change the location of BC2 to [t,2*L] and the problem is solved. I.e. both boundary conditions were defined for the same point.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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