# Trying to find values of two parameters that match the boundary conditions defined

I have been trying to solve for the values of two parameters that satisfy the boundary conditions set for a set of equations set. The below is the code.

Module[{
kbf = 0.597,
kp = 36,
ak = 7.47,
bk = 0,
ρbf = 998.2,
ρp = 3880,
cbf = 4182,
kbo = 1.38064852*10^-23,
μbf = 9.93*10^-4,
dp = 100*10^-9,
qw = 600,
cp = 773,
bulkΦ = 0.05,
Φw = 0.01,
tw = 333,
heightofchannel = 1,
pcua,
ρw,
μw,
cw,
dbw,
dtw,
ϒ,
nbt,
kw

},
{
ρw = Φw*ρp + (1 - Φw)*ρbf,
kw = kbf*(1 + ak*2*Φw + 4*Φw^2*bk),
kbf = kbf*(1 + ak*2*Φw + 4*Φw^2*bk),
μw = μbf*(1 + 2.5*Φw),
cw = (Φw*ρp*cp + (1 - Φw)*ρbf*
cbf)/ρw,
dbw = (kbo*tw)/(3*π*μbf*dp),
dtw = (0.26*kbf*μbf*Φw)/(2*kbf + kp*ρw),

ϒ = (qw*heightofchannel)/(kw*tw),
nbt = (dbw*Φw)/(dtw*ϒ),

Clear[Φ, u, t, Derivative, averagepcua],
s = ParametricNDSolve[{

k'[x] ==
1/(x^3*(2*bk*x^2 + ak*x + 1)^2)*(-2*
kbf*(-16*bk^3*x^8 - 20*ak*bk^2*x^7 - 16*bk^2*x^6 -
8*ak^2*bk*x^6 - 12*ak*bk*x^5 - ak^3*x^5 - 4*bk*x^4 -
2*ak^2*x^4 - ak*x^3)*Φ[x]),
ρ[
x] == Φ[
x]*ρp + (1 - Φ[x])* ρbf,
c[x] == (Φ[x]*ρp*
cp + (1 - Φ[x])*ρbf*cbf)/ρ[x],
μ[x] == μbf*(1 + 2.5*Φ[x]),
u'[x] == (1 - x)/(μ[x]/μw),
Φ'[x] == Φ[x]/(
nbt*(1 - (ϒ*t[x])^2))*t'[x],
t''[x] == -(
1/(k[x]/kw))*(u[x]/averagepcua +
1/kw*k'[x]*Φ'[x]*t'[x]),
Φb[x] == Integrate[u[x]*Φ[x], x]/
Integrate[u[x], x],

t == 0,
t' == 1,
u == 0,
u' == 0,
Φ == Φwp,
Φ' == 0,
k == kbf*(1 + ak*2*Φwp + 4*Φwp^2*bk)

}, {k, k', c, μ, ρ, Φ, u, t,
t', Φb}, {x, 0, 1}, {averagepcua, Φwp},
{Method -> {"Shooting"},
"StartingStepSize" -> 0.00001}

]

}]


It seems to work fine up to this point as I got the below output. However, I got error message when I try to solve for the two parameters with the below line.

FindRoot[{t'[averagepcua, Φwp] /. s,
bulkΦ - Φb[averagepcua, Φwp][
1] /. s}, {Φwp, 0.000001}]


I have inserted this line right after the ParametricNDSolve function within the Module. The below is the error message I got. I have tried to solve the issues by relating one parameter to another and reduce the number of parameters to be solved by replacing the ParametricNDSolve function with the below code.

s = ParametricNDSolve[{

k'[x] ==
1/(x^3*(2*bk*x^2 + ak*x + 1)^2)*(-2*
kbf*(-16*bk^3*x^8 - 20*ak*bk^2*x^7 - 16*bk^2*x^6 -
8*ak^2*bk*x^6 - 12*ak*bk*x^5 - ak^3*x^5 - 4*bk*x^4 -
2*ak^2*x^4 - ak*x^3)*Φ[x]),
ρ[
x] == Φ[x]*ρp + (1 - Φ[x])* ρbf,
c[x] == (Φ[x]*ρp*
cp + (1 - Φ[x])*ρbf*cbf)/ρ[x],
(*c[x]\[Equal] (Φ[x]*ρp*cp+(1-Φ[
x])*ρbf*cbf)/(Φ[x]*ρp+(1-Φ[
x])* ρbf),*)
μ[x] == μbf*(1 + 2.5*Φ[x]),
u'[x] == (1 - x)/(μ[x]/μw),
(*u'[x]\[Equal](1-x)/((μbf*(1+2.5*Φ[
x]))/μw),*)
Φ'[x] == Φ[x]/(
nbt*(1 - (ϒ*t[x])^2))*t'[x],
t''[x] == -(
1/(k[x]/kw))*(u[x]/averagepcua[x] +
1/kw*k'[x]*Φ'[x]*t'[x]),
Φb[x] == Integrate[u[x]*Φ[x], x]/
Integrate[u[x], x],
averagepcua[x] == Integrate[ρ[x]*c[x]*u[x], x],

t == 0,
t' == 1,
u == 0,
u' == 0,
Φ == Φwp,
Φ' == 0,
k == kbf*(1 + ak*2*Φwp + 4*Φwp^2*bk)

}, {k, k', c, μ, ρ, Φ, u, t,
t', Φb, averagepcua}, {x, 0, 1}, {Φwp},
{Method -> {"Shooting"},
"StartingStepSize" -> 0.00001}

]


And this for the FindRoot function

FindRoot[t'[Φwp] /. s, {Φwp, 0.01}]



None of those worked for me. I'm sorry that I don't have much knowledge about Mathematica as I'm actually quite new to it. Please do advise me on how to solve the said problem. Thanks and sorry again for the long question.

• What is the problem do you try to solve? Is it BVP for a mixture of differential and integral equations? Jan 30 at 14:25
• It is an IVP. I'm trying to solve for the value of two parameters, namely averagepcua and Φwp that gives me t'=0 at y=1. Thanks!
– Amos
Jan 30 at 14:39
• In your code Integrate[ρ[x]*c[x]*u[x], x] means NIntegrate[ρ[x1]*c[x1]*u[x1], {x1,0,x}]? Jan 30 at 16:39
• I have tried changing it to Φb[x] == NIntegrate[Integrate[u'[x1],x1]*Integrate[Φ'[x1],x1], {x1,0,x]]/ NIntegrate[Integrate[u'[x1],x1], {x1,0,x}], But I got an error message saying x1=x is not a valid limit of integration.
– Amos
Jan 31 at 3:33

In the first code we don't need to call ParametricNDSolve[] sinse we use Module. Also we don't need to compute parameters $$\rho, c, \mu$$ inside NDSolve. Therefore we can organise code as follows:

sol[fwp_, av0_] := Module[{\[CapitalPhi]wp = fwp, averagepcua = av0, kbf = 0.597, kp = 36, ak = 7.47, bk = 0, \[Rho]bf = 998.2, \[Rho]p = 3880, cbf = 4182, kbo = 1.38064852/10^23, \[Mu]bf = 9.93/10^4, dp = 100/10^9, qw = 600,
cp = 773, bulk\[CapitalPhi] = 0.05, \[CapitalPhi]w = 0.01, tw = 333, heightofchannel = 1, pcua, \[Rho]w, \[Mu]w, cw, dbw, dtw, \[CurlyCapitalUpsilon], nbt, kw}, \[Rho]w = \[CapitalPhi]w*\[Rho]p + (1 - \[CapitalPhi]w)*\[Rho]bf; kw = kbf*(1 + ak*2*\[CapitalPhi]w + 4*\[CapitalPhi]w^2*bk);
kbf = kbf*(1 + ak*2*\[CapitalPhi]w + 4*\[CapitalPhi]w^2*bk); \[Mu]w = \[Mu]bf*(1 + 2.5*\[CapitalPhi]w); cw = (\[CapitalPhi]w*\[Rho]p*cp + (1 - \[CapitalPhi]w)*\[Rho]bf*cbf)/\[Rho]w; dbw = (kbo*tw)/(3*Pi*\[Mu]bf*dp); dtw = (0.26*kbf*\[Mu]bf*\[CapitalPhi]w)/(2*kbf + kp*\[Rho]w);
\[CurlyCapitalUpsilon] = (qw*heightofchannel)/(kw*tw); nbt = (dbw*\[CapitalPhi]w)/(dtw*\[CurlyCapitalUpsilon]); \[Rho] = \[CapitalPhi][x]*\[Rho]p + (1 - \[CapitalPhi][x])*\[Rho]bf; c = (\[CapitalPhi][x]*\[Rho]p*cp + (1 - \[CapitalPhi][x])*\[Rho]bf*cbf)/\[Rho]; \[Mu] = \[Mu]bf*(1 + 2.5*\[CapitalPhi][x]); x0 = $MachineEpsilon; s = NDSolve[{Derivative[k][x] == (1/(x^3*(2*bk*x^2 + ak*x + 1)^2))*(-2*kbf*(-16*bk^3*x^8 - 20*ak*bk^2*x^7 - 16*bk^2*x^6 - 8*ak^2*bk*x^6 - 12*ak*bk*x^5 - ak^3*x^5 - 4*bk*x^4 - 2*ak^2*x^4 - ak*x^3)*\[CapitalPhi][x]), Derivative[u][x] == (1 - x)/(\[Mu]/\[Mu]w), Derivative[\[CapitalPhi]][x] == (\[CapitalPhi][x]/(nbt*(1 - (\[CurlyCapitalUpsilon]*t[x])^2)))*Derivative[t][x], Derivative[t][x] == (-(1/(k[x]/kw)))*(u[x]/averagepcua + (1/kw)*Derivative[k][x]*Derivative[\[CapitalPhi]][x]*Derivative[t][x]), t[x0] == 0, Derivative[t][x0] == 1, u[x0] == 0, \[CapitalPhi][x0] == \[CapitalPhi]wp, k[x0] == kbf*(1 + ak*2*\[CapitalPhi]wp + 4*\[CapitalPhi]wp^2*bk)}, {k, \[CapitalPhi], u, t}, {x, x0, 1}]; s]  Now we can evaluate, for example, sol[.01, .01] and plot u, k' as {Plot[u[x] /. s, {x, x0, 1}, AxesLabel -> {"x", "u"}], Plot[t'[x] /. s, {x, x0, 1}, AxesLabel -> {"x", "k'"}]} Practically we get with this parameters desired solution since t'[1.] /. s Out[]= {-0.000362355}  Also some numerical integration can be implemented as follows \[CapitalPhi]b[x1_?NumericQ] := NIntegrate[u[x]*\[CapitalPhi][x] /. s, {x, x0, x1}]/ NIntegrate[u[x] /. s, {x, x0, x1}] \[CapitalPhi]b[.5] {3.345431280729773}  With Module[] we can also compute unknown parameters \[CapitalPhi]wp, averagepcua by using additional conditions u'==0,t'==0. For this we define a new function sol1[fwp_, av0_] := Module[{\[CapitalPhi]wp = fwp, averagepcua = av0, kbf = 0.597, kp = 36, ak = 7.47, bk = 0, \[Rho]bf = 998.2, \[Rho]p = 3880, cbf = 4182, kbo = 1.38064852/10^23, \[Mu]bf = 9.93/10^4, dp = 100/10^9, qw = 600, cp = 773, bulk\[CapitalPhi] = 0.05, \[CapitalPhi]w = 0.01, tw = 333, heightofchannel = 1, pcua, \[Rho]w, \[Mu]w, cw, dbw, dtw, \[CurlyCapitalUpsilon], nbt, kw,u1,t1}, \[Rho]w = \[CapitalPhi]w*\[Rho]p + (1 - \[CapitalPhi]w)*\[Rho]bf; kw = kbf*(1 + ak*2*\[CapitalPhi]w + 4*\[CapitalPhi]w^2*bk); kbf = kbf*(1 + ak*2*\[CapitalPhi]w + 4*\[CapitalPhi]w^2*bk); \[Mu]w = \[Mu]bf*(1 + 2.5*\[CapitalPhi]w); cw = (\[CapitalPhi]w*\[Rho]p*cp + (1 - \[CapitalPhi]w)*\[Rho]bf*cbf)/\[Rho]w; dbw = (kbo*tw)/(3*Pi*\[Mu]bf*dp); dtw = (0.26*kbf*\[Mu]bf*\[CapitalPhi]w)/ (2*kbf + kp*\[Rho]w); \[CurlyCapitalUpsilon] = (qw*heightofchannel)/(kw*tw); nbt = (dbw*\[CapitalPhi]w)/(dtw*\[CurlyCapitalUpsilon]); \[Rho] = \[CapitalPhi][x]*\[Rho]p + (1 - \[CapitalPhi][x])*\[Rho]bf; c = (\[CapitalPhi][x]*\[Rho]p*cp + (1 - \[CapitalPhi][x])*\[Rho]bf*cbf)/\[Rho]; \[Mu] = \[Mu]bf*(1 + 2.5*\[CapitalPhi][x]); x0 =$MachineEpsilon;
{u1,t1}={u',t'}/. NDSolve[{Derivative[k][x] ==
-((1/(x^3*(2*bk*x^2 + ak*x + 1)^2))*(2*kbf*(-16*bk^3*x^8 -
20*ak*bk^2*x^7 - 16*bk^2*x^6 - 8*ak^2*bk*x^6 -
12*ak*bk*x^5 - ak^3*x^5 - 4*bk*x^4 - 2*ak^2*x^4 - ak*x^3)*
\[CapitalPhi][x])), Derivative[u][x] == (1 - x)/(\[Mu]/\[Mu]w),
Derivative[\[CapitalPhi]][x] == (\[CapitalPhi][x]*Derivative[t][x])/
(nbt*(1 - (\[CurlyCapitalUpsilon]*t[x])^2)), Derivative[t][x] ==
-((u[x]/averagepcua + (Derivative[k][x]*Derivative[\[CapitalPhi]][x]*
Derivative[t][x])/kw)/(k[x]/kw)), t[x0] == 0,
Derivative[t][x0] == 1, u[x0] == 0, \[CapitalPhi][x0] == \[CapitalPhi]wp,
k[x0] == kbf*(1 + ak*2*\[CapitalPhi]wp + 4*\[CapitalPhi]wp^2*bk)}, {k, \[CapitalPhi], u, t},
{x, x0, 1}][]; {u1,t1}]


This function can be evaluated as follows

ff[p1_?NumericQ, p2_?NumericQ] := sol1[p1, p2]

sf=FindRoot[ff[p1, p2], {{p1, .01}, {p2, .01}}]


FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

Out[]= {p1 -> 0.0260872, p2 -> 0.0206202}


This is the best solution what we can get with NDSolve[]

ff[p1, p2] /.sf

Out[]= {3.67241*10^-11, -0.0000885718}


It looks as it shown above, and unknown parameters are {\[CapitalPhi]wp, averagepcua}={0.0260872, 0.0206202}

Finally we can also use ParametricNDSolve[] as follows

par = {kbf = 0.597, kp = 36, ak = 7.47, bk = 0, \[Rho]bf = 998.2, \[Rho]p = 3880, cbf = 4182, kbo = 1.38064852/10^23,
\[Mu]bf = 9.93/10^4, dp = 100/10^9, qw = 600, cp = 773, bulk\[CapitalPhi] = 0.05, \[CapitalPhi]w = 0.01, tw = 333, heightofchannel = 1};
\[Rho]w = \[CapitalPhi]w*\[Rho]p + (1 - \[CapitalPhi]w)*\[Rho]bf; kw = kbf*(1 + ak*2*\[CapitalPhi]w + 4*\[CapitalPhi]w^2*bk);
kbf = kbf*(1 + ak*2*\[CapitalPhi]w + 4*\[CapitalPhi]w^2*bk); \[Mu]w = \[Mu]bf*(1 + 2.5*\[CapitalPhi]w);
cw = (\[CapitalPhi]w*\[Rho]p*cp + (1 - \[CapitalPhi]w)*\[Rho]bf*cbf)/\[Rho]w;
dbw = (kbo*tw)/(3*Pi*\[Mu]bf*dp); dtw = (0.26*kbf*\[Mu]bf*\[CapitalPhi]w)/(2*kbf + kp*\[Rho]w); \[CurlyCapitalUpsilon] = (qw*heightofchannel)/(kw*tw);
nbt = (dbw*\[CapitalPhi]w)/(dtw*\[CurlyCapitalUpsilon]); \[Rho] = \[CapitalPhi][x]*\[Rho]p + (1 - \[CapitalPhi][x])*\[Rho]bf;
c = (\[CapitalPhi][x]*\[Rho]p*cp + (1 - \[CapitalPhi][x])*\[Rho]bf*cbf)/\[Rho];
\[Mu] = \[Mu]bf*(1 + 2.5*\[CapitalPhi][x]); x0 = \$MachineEpsilon;
eq = {Derivative[k][x] == -((2*kbf*(-16*bk^3*x^8 - 20*ak*bk^2*x^7 - 16*bk^2*x^6 - 8*ak^2*bk*x^6 - 12*ak*bk*x^5 -
ak^3*x^5 - 4*bk*x^4 - 2*ak^2*x^4 - ak*x^3)*\[CapitalPhi][x])/(x^3*(2*bk*x^2 + ak*x + 1)^2)),
Derivative[u][x] == (1 - x)/(\[Mu]/\[Mu]w), Derivative[\[CapitalPhi]][x] == (\[CapitalPhi][x]*Derivative[t][x])/
(nbt*(1 - (\[CurlyCapitalUpsilon]*t[x])^2)), Derivative[t][x] ==
-((u[x]/averagepcua + (Derivative[k][x]*Derivative[\[CapitalPhi]][x]*Derivative[t][x])/kw)/(k[x]/kw)), t[x0] == 0,
Derivative[t][x0] == 1, u[x0] == 0, \[CapitalPhi][x0] == \[CapitalPhi]wp, k[x0] == kbf*(1 + ak*2*\[CapitalPhi]wp + 4*\[CapitalPhi]wp^2*bk)};


Solve and evaluate unknown parameters

sol1 =
ParametricNDSolveValue[
eq, {u', t'}, {x, x0, 1}, {\[CapitalPhi]wp, averagepcua}];

ff[p1_?NumericQ, p2_?NumericQ] := sol1[p1, p2]

sf = FindRoot[ff[p1, p2], {{p1, .026}, {p2, .02}}] // Quiet

(*Out[]= {p1 -> 0.0252197, p2 -> 0.0199972}*)

{\[CapitalPhi]wp, averagepcua} = {p1, p2} /. sf;


Plot numerical solution

s = NDSolve[eq, {k, \[CapitalPhi], u, t}, {x, x0, 1}]

{Plot[u'[x] /. s, {x, x0, 1}, AxesLabel -> {"x", "u'"}],
Plot[t'[x] /. s, {x, x0, 1}, AxesLabel -> {"x", "t'"}]}

{Plot[u[x] /. s, {x, x0, 1}, AxesLabel -> {"x", "u"}],
Plot[t[x] /. s, {x, x0, 1}, AxesLabel -> {"x", "t"}]}
` • Hello. Thanks for taking time to answer my question. I have read through the coding. Actually the reason why I used ParametricNDSolve is because the value of Φwp and averagepcua are not known initially. I'm actually trying to find the right combination of Φwp and averagepcua that give t'=0 at x=1. Please correct me if I'm wrong but I think NDSolve only helps in cases where the value of all parameters are known and only numerical analysis is required.
– Amos
Feb 1 at 14:52
• By the way, I'm currently trying to solve for the two parameters with ParametricNDSolve and a While loop. I have removed all the terms without derivatives and substituted the said terms with their equations into the equations such as t'' where they are used in the ParametricNDSolve function. Then I use a While loop to keep running the ParametricNDSolve with different value of Φwp and averagepcua until I get what I want. I'm not sure whether if it works yet as it is still running now. Please do let me know if you think this will works and a little suggestions thanks!
– Amos
Feb 1 at 14:57
• @Amos See update to my answer. Feb 1 at 18:15
• Hello, sorry for replying late, have been busy the past few days. I'm currently trying to solve it in the way you did. I'm facing a few errors here which I think it's due to the way I wrote the coding. I will try to solve the errors first. Thanks!
– Amos
Feb 7 at 12:43
• It works, thank you very much!
– Amos
Feb 9 at 5:25