I am trying to solve a multiphysic problem using finite elements and a Newton Raphson solution scheme. I have two non-linear subsystems that are coupled bi-directionally.

The first subsystem includes an equation that I solve for $p_f$:

$\nabla u - \nabla a^3 \nabla p_f = \nabla a^3 B$

This equation becomes very non-linear, due to the parameter $a$ with an exponent of 3! $B$ are body forces and $u$ is the solution of the second subsystem that I am trying to solve. This second subsystem includes the Stokes equation that depends on the gradient of $p_f$. Therefore, both system are coupled bi-directionally.

I am using a Newton-Raphson solution scheme for each subsystem. First, I solve for $p_f$. After convergence occurred, I pass this solution to the second subsystem, which I solve using a Newton Raphson scheme, too. Then, the solution $u$ is passed to the first subsystem ... until the solutions of both subsystems do not change anymore. See the flow chart below.

Flow Chart

This procedure works very well so far. However, if the values of the parameter $a$ become very low (thus the gradient of $a$ becomes very, very low), there is no convergence anymore! The reason seems to be that the coupling between both subsystem becomes very strong. A small change in one subsystem creates a significant changes in the other one (especially where the gradient of $a$ is very low). I tried a variety of good initial conditions, however, since the first subsystem reacts very sensitiv, convergence has been impossible so far.

In the figure below, you can see how the solution of this approach ("Segregated Approach"-blue line) is different from the "true" solution ("Fully Coupled"-red line). The "true solution" is the result of a benchmark, which I calculated by solving a single system of equations that includes the equations of both subsystems.

Differences in solution

The greatest differences between both solutions are located where the gradient of $a$ becomes very small. See the figure below (in this figure the values of $p_f$ and $a$ are normalized to range between 0 and 1):

Parameter $a$ and the 'wrong' solution for $p_f$

So: Is a Newton-Raphson solution scheme the wrong method for solving a strongly coupled system, if there are low gradients in a subsystem? Could my method be fixed to solve the problem that I described? If yes, how? Are there alternative solution schemes for this problem?

  • $\begingroup$ It's probably the segregation scheme that's killing you not the inner non-linear solver. Why don't you want to solve the fully coupled problem? $\endgroup$
    – Bill Barth
    Commented Mar 5, 2015 at 14:18
  • $\begingroup$ 1. The model that I'm modifying makes it very hard to add a new subsystem to the existing system of equations 2. Later, the model should run on a larger scale, and it seemed less time and memory consuming to use a segregation scheme $\endgroup$
    – Johann
    Commented Mar 5, 2015 at 14:23
  • $\begingroup$ Find a new program to extend? I think the odds are good that the segregated version is converging to the wrong answer. If the inner Newton steps appear converged, and the outer iterations between the segregated models appear converged, then segregation may not be working. How many different initial guesses have you tried? It's possible that the fully coupled version has an easier time hitting what you think of as the right solution. How did you generate the fully coupled solution if the model is hard to extend, btw? Have you tried the method of manufactured solutions on both models? $\endgroup$
    – Bill Barth
    Commented Mar 5, 2015 at 16:27
  • $\begingroup$ The fully coupled solution is from a simple model I wrote in MATLAB (much too simple for later usage, but good enough to run this benchmark). I used the solution of the fully-coupled code as initial condition and did not get any convergence! I have not tried the method of manufactured solutions yet. $\endgroup$
    – Johann
    Commented Mar 5, 2015 at 22:06
  • 1
    $\begingroup$ The nonlinear term is $\nabla (a^3)$ or $(\nabla a)^3$? $\endgroup$
    – nicoguaro
    Commented Mar 19, 2015 at 22:26

3 Answers 3


The issues you're running into now are not a failing of Newton-Raphson, but a question of coupling. You're doing iterated sequential coupling -- solving each equation sequentially and then iterating until (hopeful) convergence.

No solver choice in place of NR is going to fix this lack of convergence, as long as you are doing iterated sequential coupling.

Instead, you want to consider different coupling schemes that use more information about the coupling terms. What @paul suggests is effectively equivalent to doing a quasi-Newton method using an approximation of the Jacobian of the coupled system that is block diagonal. Note that this is NOT Newton-Raphson any more, because you do not have the true Jacobian of the coupled system -- instead you form an approximation of the Jacobian using fewer terms than the true Jacobian and use this instead. This is a good starting point because it requires no additional knowledge of the physics -- you can build this entirely from pieces you already have. Norms for the coupled system become a combination of the norms for the two systems you have. Residuals for the coupled system become the vector of both systems' residuals. And the approximate inverse is block diagonal.

It is possible, but unlikely, that this will help much. The issue is that the approximate (block diagonal) inverse will be too far from the true inverse. One potential option is to try Jacobian-Free Newton Krylov to get a better inverse. This is likely to not work all that much better until you start to add off-diagonal blocks to the approximate inverse.

From a practical standpoint -- your best bet is to bite the bullet and find a way to start to include at least some information about how the equations couple in your approximate inverse, i.e. including off-diagonal blocks in your approximate Jacobian for use in a fully implicit coupling scheme.


Instead of running Newton-Raphson to convergence of each subsystem, try 1 iteration on the first subsystem followed by 1 iteration on the second subsystem. This may keep the subsystems more coupled, not going to "distant" unrelated sub solutions. Repeat this two step iteration and see if it converges.


This is a short answer and mostly a workaround. This could hit you back again later on, but can provide you with a temporary solution. You can try to relax the solution to your Newton-Raphson solver. What I mean by this is do not take the full solution of your non-linear system, but just take a factor of the previous solution and of the new solution (say 0.1 new solution, 0.9 old solution), and use this intermediate solution in your other solver. This means that you would be relaxing your global system of equation (and is similar conceptually to the above answer of using a single Newton-Raphson iteration).

This usually leads to more robust convergence, but it can turn out to be extremely slow.


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.