fronts.solve_from_guess(D, i, b, o_guess, guess, radial=False, max_nodes=1000, verbose=0)

Alternative solver for problems with a Dirichlet boundary condition.

Given a positive function D, scalars \(\theta_i\), \(\theta_b\) and \(o_b\), and coordinate unit vector \(\mathbf{\hat{r}}\), finds a function \(\theta\) of r and t such that:

\[\begin{split}\begin{cases} \dfrac{\partial\theta}{\partial t} = \nabla\cdot\left[D(\theta))\dfrac{\partial\theta}{\partial r} \mathbf{\hat{r}}\right ] & r>r_b(t),t>0\\ \theta(r, 0) = \theta_i & r>0 \\ \theta(r_b(t), t) = \theta_b & t>0 \\ r_b(t) = o_b\sqrt t \end{cases}\end{split}\]

Alternative to the main solve() function. This function requires a starting mesh and guess of the solution. It is significantly less robust than solve(), and will fail to converge in many cases that the latter can easily handle (whether it converges will usually depend heavily on the problem, the starting mesh and the guess of the solution; it will raise a RuntimeError on failure). However, when it converges it is usually faster than solve(), which may be an advantage for some use cases. You should nonetheless prefer solve() unless you have a particular use case for which you have found this function to be better.

Possible use cases include refining a solution (note that solve() can do that too), optimization runs in which known solutions make good first approximations of solutions with similar parameters and every second of computing time counts, and in the implementation of other solving algorithms. In all these cases, solve() should probably be used as a fallback for when this function fails.

  • D (callable or sympy.Expression or str or float) –

    Callable that evaluates \(D\) and its derivatives, obtained from the fronts.D module or defined in the same manner—i.e.:

    • D(theta) evaluates and returns \(D\) at theta

    • D(theta, 1) returns both the value of \(D\) and its first derivative at theta

    • D(theta, 2) returns the value of \(D\), its first derivative, and its second derivative at theta

    where theta may be a single float or a NumPy array.

    Alternatively, instead of a callable, the argument can be the expression of \(D\) in the form of a string or sympy.Expression with a single variable. In this case, the solver will differentiate and evaluate the expression as necessary.

  • i (float) – Initial condition, \(\theta_i\).

  • b (float) – Imposed boundary value, \(\theta_b\).

  • o_guess (numpy.array_like, shape (n_guess,)) –

    Starting mesh in terms of the Boltzmann variable. Must be strictly increasing. o_guess[0] is taken as the value of the parameter \(o_b\), which determines the behavior of the boundary. If zero, it implies that the boundary always exists at \(r=0\). It must be strictly positive if radial is not False. A non-zero value implies a moving boundary.

    On the other end, o_guess[-1] must be large enough to contain the solution to the semi-infinite problem.

  • guess (float or numpy.array_like, shape (n_guess,)) – Starting guess of the solution at the points in o_guess. If a single value, the guess is assumed uniform.

  • radial ({False, 'cylindrical', 'polar', 'spherical'}, optional) –

    Choice of coordinate unit vector \(\mathbf{\hat{r}}\). Must be one of the following:

    • False (default): \(\mathbf{\hat{r}}\) is any coordinate unit vector in rectangular (Cartesian) coordinates, or an axial unit vector in a cylindrical coordinate system

    • 'cylindrical' or 'polar': \(\mathbf{\hat{r}}\) is the radial unit vector in a cylindrical or polar coordinate system

    • 'spherical': \(\mathbf{\hat{r}}\) is the radial unit vector in a spherical coordinate system

  • max_nodes (int, optional) – Maximum allowed number of mesh nodes.

  • verbose ({0, 1, 2}, optional) –

    Level of algorithm’s verbosity. Must be one of the following:

    • 0 (default): work silently

    • 1: display a termination report

    • 2: also display progress during iterations


solution – See Solution for a description of the solution object. Additional fields specific to this solver are included in the object:

  • o (numpy.ndarray, shape (n,)) – Final solver mesh, in terms of the Boltzmann variable.

  • niter (int) – Number of iterations required to find the solution.

  • rms_residuals (numpy.ndarray, shape (n-1,)) – RMS values of the relative residuals over each mesh interval.

Return type:


See also



This function works by transforming the partial differential equation with the Boltzmann transformation using ode() and then solving the resulting ODE with SciPy’s collocation-based boundary value problem solver scipy.integrate.solve_bvp() and a two-point Dirichlet condition that matches the boundary and initial conditions of the problem. Upon that solver’s convergence, it runs a final check on whether the candidate solution also satisfies the semi-infinite condition (which implies \(d\theta/do\to0\) as \(o\to\infty\)).