Solving Nonlinear Problems
FlexPDE automatically recognizes when a problem is nonlinear and modifies its strategy accordingly. The solution method used by FlexPDE is a modified Newton-Raphson iteration procedure. This is a "descent" method, which tries to fall down the gradient of an energy functional until minimum energy is achieved (i.e. the gradient of the functional goes to zero). If the functional is nearly quadratic, as it is in simple diffusion problems, then the method converges quadratically (the relative error is squared on each iteration). The default strategy implemented in FlexPDE is frequently sufficient to determine a solution without user intervention. But in cases of strong nonlinearities, it may be necessary for the user to help guide FlexPDE to a valid solution. There are several techniques that can be used to help the solution process.
In nonlinear time-dependent problems, the default behavior is to take a single Newton step at each timestep, on the assumption that any nonlinearities will be sensed by the timestep controller, and that timestep adjustments will guarantee an accurate evolution of the system from the given initial conditions. In this mode, the derivatives of the solution with respect to the variables is computed once at the beginning of the timestep, and are not updated.
In the case of nonlinear steady-state problems, the situation is somewhat more complicated. We are not guaranteed that the system will have a unique solution, and even if it does, we are not guaranteed that FlexPDE will be able to find it.
|•||Start with a Good Initial Value|
Providing an initial value which is near the correct solution will aid enormously in finding a solution. Be particularly careful that the initial value matches the boundary conditions. If it does not, serious excursions may be excited in the trial solution, leading to solution difficulties.
|•||Use STAGES to Gradually Activate the Nonlinear Terms|
You can use the staging facility of FlexPDE to gradually increase the strength of the nonlinear terms. Start with a linear or nearly linear system, and allow FlexPDE to find a solution which is consistent with the boundary conditions. Then use this solution as a starting point for a more strongly nonlinear system. By judicious use of staging, you can creep up on a solution to very nasty problems.
|•||Use artificial diffusion to stabilize solutions|
Gibbs phenomena are observed in signal processing when a discontinuous signal is reconstructed from its Fourier components. The characteristic of this phenomenon is ringing, with overshoots and undershoots in the recovered signal. Similar phenomena can be observed in finite element models when a sharp transition is modeled with an insufficient density of mesh cells. In signal processing, the signal can be smoothed by use of a "window function". This is essentially a low-pass filter that removes the high frequency components of the signal. In partial differential equations, the diffusion operator Div(grad(u)) is a low-pass filter that can be used to smooth oscillations in the solution. See the Technical Note "Smoothing Operators in PDE's" for technical discussion of this operator. In brief, you can use a term eps*Div(Grad(u)) in a PDE to smooth oscillations of spatial extent D by setting eps=D^2/pi^2 in steady state or eps=2*D*c/pi in time dependence (where c is the signal propagation velocity). The term should also be scaled as necessary to provide dimensional consistency with the rest of the terms of the equation. Use of such a term merely limits the spatial frequency components of the solution to those which can be adequately resolved by the finite element mesh.
|•||Use CHANGELIM to Control Modifications|
The selector CHANGELIM limits the amount by which any nodal value in a problem may be modified on each Newton-Raphson step. As in a one-dimensional Newton iteration, if the trial solution is near a local maximum of the functional, then shooting down the gradient will try to step an enormous distance to the next trial solution. FlexPDE applies a backtracking algorithm to try to find the step size of optimal residual reductions, but it also limits the size of each nodal change to be less than CHANGELIM times the average value of the variable. The default value for CHANGELIM is 0.5, but if the initial value (or any intermediate trial solution) is sufficiently far from the true solution, this value may allow wild excursions from which FlexPDE is unable to recover. Try cutting CHANGELIM to 0.1, or in severe cases even 0.01, to force FlexPDE to creep toward a valid solution. In combination with a reasonable initial value, even CHANGELIM=0.01 can converge in a surprisingly short time. Since CHANGELIM multiplies the RMS average value, not each local value, its effect disappears when a solution is reached, and quadratic final convergence is still achieved.
|•||Watch Out for Negative Values|
FlexPDE uses piecewise polynomials to approximate the solution. In cases of rapid variation of the solution over a single cell, you will almost certainly see severe under-shoot in early stages. If you are assuming that the value of your variable will remain positive, don't. If your equations lose validity in the presence of negative values, perhaps you should recast the equations in terms of the logarithm of the variable. In this case, even though the logarithm may go negative, the implied value of your actual variable will remain positive.
|•||Recast the Problem in a Time-Dependent Form|
Any steady-state problem can be viewed as the infinite-time limit of a time-dependent problem. Rewrite your PDE's to have a time derivative term which will push the value in the direction of decreasing deviation from solution of the steady-state PDE. (A good model to follow is the time-dependent diffusion equation DIV(K*GRAD(U)) = DT(U). A negative value of the divergence indicates a local maximum in the solution, and results in driving the value downward.) In this case, "time" is a fictitious variable analogous to the "iteration count" in the steady-state N-R iteration, but the time-dependent formulation allows the timestep controller to guide the evolution of the solution.