1
$\begingroup$

As a very beginner into Mathematica, I encounter error message as shown below when I am solving the system of three second order ODEs numerically using ParametricNDSolveValue:

FindRoot::dfmin: The minimal damping factor of 1/10000 has been reached.

ParametricNDSolveValue::fempsf: PDESolve could not find a solution.

Below is the essential part of my codes, starting with list of constants,

(*List of constants*)
w = 1; 
p = 1; 
q = 0.05; 
u = 1; 
Rba = 0; 
Rca = 1; 
Kb = 1.38*10^-23; 
T = 25 + 273; 
miumaxC = 0;
alphaA = 20;
alphaB = 20;
alphaC = 1;
betaA = 1;
betaB = 1;
betaC = 1;
kgA = 1;
kgB = 1;
kgC = 1;

Intermediate functions,

(*intermediate functions*)
cp1[x_] = miumaxA/2*(1 - Tanh[alphaA (x - betaA)])*Kb*T; 
cp2[x_] = miumaxB/2*(1 - Tanh[alphaB (x - betaB)])*Kb*T; 
cp3[x_] = miumaxC/2*(1 - Tanh[alphaC (x - betaC)])*Kb*T; 
diffintermediate1[x_] = a[x]*x^2*cp1'[x];
diffintermediate2[x_] = b[x]*x^2*cp2'[x];
diffintermediate3[x_] = (phi[x] - 1)*x^2*cp3'[x];

Defining ODEs and boundary conditions before solving,

(*Dimensionless species balance equations with neumann BCs for \
species A,B,and C*)
ODE1 = a''[x] + 2/x*a'[x] + 1/(Kb*T)/x^2*diffintermediate1'[x] - 
   TM^2*phi[x]*a[x] == 
  NeumannValue[0, x == 0] + 
   NeumannValue[kgA*(a[x] - 1), x == 1]; ODE2  = 
 b''[x] + 2/x*b'[x] + 1/(Kb*T)/x^2*diffintermediate2'[x] + 
   TM^2*phi[x]/p*(a[x] - w*b[x]) == 
  NeumannValue[0, x == 0] + NeumannValue[kgB/p* (b[x]), x == 1];
ODE3  = q*(phi''[x] + 2/x*phi'[x] + 
       1/(Kb*T)/x^2*diffintermediate3'[x]) - TM^2*u*w*phi[x]*b[x] == 
   NeumannValue[0, x == 0] + NeumannValue[kgC/q* (phi[x] - 1), x == 1];

(*Grouping the ODEs for input in NDSolveValue function*)
ODEs = {ODE1, ODE2 , ODE3}; 

(*Solving the system of ODES*)
solutions = 
  ParametricNDSolveValue[{ODEs}, {a, b, phi}, {x, 0, 1}, {TM, 
    miumaxA, miumaxB}]; 

I would get the abovementioned error message when I try to calculate and plot the solutions at TMplot = 10 to investigate the effect of different miumax. Below is the code for plotting the solutions:

MiumaxList = {-4,-2,0,2,4};
(*Solution profile plotting at different miumaxA*)
TMplot  =10;

(*Creating legend colours*)
legendColours = {};
legendEntries = {};
For[i=1,i<=Length[MiumaxList],i++,
legendColours = Insert[legendColours,RGBColor[0.2*i,0.35,0.7],i];
]

(*Creating list of solutions for varying Miumax*)
solutionListA = {};
solutionListB = {};
solutionListPhi = {};
For[i=1,i<=Length[MiumaxList],i++,
solutionListA= Insert[solutionListA,solutions[TMplot,MiumaxList[[i]],MiumaxList[[i]]][[1]][x] ,i];
solutionListB = Insert[solutionListB,solutions[TMplot,MiumaxList[[i]],MiumaxList[[i]]][[2]][x] ,i];
solutionListPhi = Insert[solutionListPhi,solutions[TMplot,MiumaxList[[i]],MiumaxList[[i]]][[3]][x] ,i];
];

(*Plotting of solutions a, b, and phi for varying Miumax*)
Plot[solutionListA,{x,0,1},PlotRange->{Automatic,Full},PlotLabel->"a (x) at different miumax pair, TM = "[TMplot],AxesLabel->{"x","a(x)"},PlotStyle->legendColours,PlotLegends->LineLegend[legendColours,MiumaxList]]
Plot[solutionListB,{x,0,1},PlotRange->{Automatic,Full},PlotLabel->"b (x) at different miumax pair, TM = "[TMplot],AxesLabel->{"x","a(x)"},PlotStyle->legendColours,PlotLegends->LineLegend[legendColours,MiumaxList]]
Plot[solutionListPhi,{x,0,1},PlotRange->{Automatic,Full},PlotLabel->"a (x) at different miumax pair, TM = "[TMplot],AxesLabel->{"x","a(x)"},PlotStyle->legendColours,PlotLegends->LineLegend[legendColours,MiumaxList]]
$\endgroup$

1 Answer 1

2
$\begingroup$

There are two things you can try. In version 13.1 parametric functions when used with the finite element method can take a new initial seeding. This is relevant in your case as the nonliterary is very strong as numbers increase and a good initial seed might help to push the solutions further out.

Let's have a look:

(*Solving the system of ODES*)
solutions = 
  ParametricNDSolveValue[{ODEs}, {a, b, phi}, {x, 0, 1}, {TM, miumaxA,
     miumaxB}];

We can compute a solution for this set of parameters:

oldSolution = solutions[3, -4, -4]
(* InterpolatingFunction[.... *)

But when we increase the 3 to a 4 we get a convergence error:

solutions[4, -4, -4]

enter image description here

What you can do now is to use the oldSolution as an initial seed for the next solution like so:

solutions[4, -4, -4, 
 "InitialSeeding" -> 
  Thread[Through[{a, b, phi}[x]] == Through[oldSolution[x]]]]

I can not say how far you can push these out but give it a shot. If this does not work, not all hope is lost. There is a second way to do this. In that case you convert your stationary ODEs to time dependent ODEs and time integrate them for a pseudo time from 0 to 1 and use for example the first parameter like

NDSolve[ tsolutions[t*10, -4, -4],..., {t,0,1},...]

This has the advantage that the adaptive time stepping algorithm of NDSolve will take care of the step size selection; on the downside you need to set up a time dependent problem and initial conditions.

Both approaches are shown in the monograph on Hyperelasticity.

Good luck!

$\endgroup$
4
  • $\begingroup$ Thank you so much for the advice. I just tried the first method, which is to alter the initial seeding by calling it using solutions[4, -4, -4, "InitialSeeding" -> Thread[Through[{a, b, phi}[x]] == Through[oldSolution[x]]]]. However, I would get error: "ParametricNDSolveValue::fpct: Too many parameters in {TM,miumaxA,miumaxB} to be filled from ... ..." It seems like I have to change the initial seedings in ParametricNDSolveValue[... ...]. Is there any workaround that changes the initial seeding after declaring ParametricNDSolveValue? $\endgroup$
    – Johnson
    Commented Jul 11, 2022 at 9:48
  • $\begingroup$ @joe, what is you $Version? $\endgroup$
    – user21
    Commented Jul 11, 2022 at 10:48
  • $\begingroup$ I am currently using ver 13.0.1.0 $\endgroup$
    – Johnson
    Commented Jul 11, 2022 at 12:22
  • 1
    $\begingroup$ @joe, as mentioned in the beginning of the post, the above works in version 13.1. In version 13.0 you'd have to try the second method that I outlined. Make this a time dependent problem and add a factor of t in front of the parameters. $\endgroup$
    – user21
    Commented Jul 11, 2022 at 13:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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