# Solving non linear system of equations with Newton Raphson method
**URL:** https://community.freefem.org/t/solving-non-linear-system-of-equations-with-newton-raphson-method/3733
**Category:** General Discussion
**Created:** [February 6, 2025, 5:41pm UTC](https://community.freefem.org/t/solving-non-linear-system-of-equations-with-newton-raphson-method/3733 "2025-02-06T17:41:10Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@Try\_everything](https://community.freefem.org/u/Try_everything)
#### Post date: [February 6, 2025, 5:41pm UTC](https://community.freefem.org/t/solving-non-linear-system-of-equations-with-newton-raphson-method/3733/1 "2025-02-06T17:41:10Z")
Hi , I am trying to solve steady boussinesq approximation of governing equations of flow compromising conservation of mass, momentum, energy. I weak formulated the equation and I need to use Newton raphson method to solve the non-linear residuals. I am requesting some guidance in how to solve problem of multiple coupled non linear residuals.
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [February 7, 2025, 9:15am UTC](https://community.freefem.org/t/solving-non-linear-system-of-equations-with-newton-raphson-method/3733/2 "2025-02-07T09:15:46Z")
You have an example with comments in the FreeFem documentation p.75 :
Newton Method for the Steady Navier-Stokes equations
> **[FreeFEM-documentation.pdf](https://doc.freefem.org/pdf/FreeFEM-documentation.pdf)**
>
> 33.85 MB
---
### Author:  [@Try\_everything](https://community.freefem.org/u/Try_everything)
#### Post date: [February 13, 2025, 2:45pm UTC](https://community.freefem.org/t/solving-non-linear-system-of-equations-with-newton-raphson-method/3733/3 "2025-02-13T14:45:16Z")
Thanks for your reply. I tried it. But I am implemented the following code and I made an error here. How to correct and identify it?
```auto
// Parameters
real Pr = 0.7;
real Ra = 1000;
real gamma = 100; // penalty parameter
real eps = 1e-5;
// Mesh
mesh Th = square(33, 33); // Labels: left=1, right=2, bottom=3, top=4
plot(Th, wait=true);
// Fespace
fespace Xh(Th, P2);
Xh u1, u2, v1, v2, du1, du2, u1T, u2T;
fespace Mh(Th, P1);
Mh T, q, dT, pp;
// Stream Function Space
fespace Ph(Th, P2);
Ph psi, phi;
// Macros
macro Grad(u1, u2) [dx(u1), dy(u1), dx(u2), dy(u2)] //
macro UgradV(u1, u2, v1, v2) [[u1, u2]'*[dx(v1), dy(v1)], [u1, u2]'*[dx(v2), dy(v2)]] //
macro div(u1, u2) (dx(u1) + dy(u2)) //
macro DivInner(u1, u2, v1, v2) ((dx(u1) + dy(u2)) * (dx(v1) + dy(v2))) //
macro Graddotgrad(T, q) (dx(q) * dx(T) + dy(q) * dy(T)) //
macro Advect(u1, u2, T) (u1 * dx(T) + u2 * dy(T)) //
// Initialization
u1 = 0;
u2 = 0;
T = 0;
// Viscosity loop
while (1) {
int n;
real err = 0;
// Newton loop
for (n = 0; n < 15; n++) {
// Solve velocity equation
solve velocity ([du1, du2], [v1, v2])
= int2d(Th)(
Pr * (Grad(du1, du2)' * Grad(v1, v2))
+ UgradV(du1, du2, u1, u2)' * [v1, v2]
+ UgradV(u1, u2, du1, du2)' * [v1, v2]
+ gamma * DivInner(du1, du2, v1, v2) // divergence penalty
- Ra * Pr * dT * v2
)
-int2d(Th)( Pr * (Grad(u1, u2)' * Grad(v1, v2))
+ UgradV(u1, u2, u1, u2)' * [v1, v2]
+ gamma * DivInner(u1, u2, v1, v2) // divergence penalty
- Ra * Pr * T * v2
)
+ on(1, u1=0, u2=0)
+ on(2, u1=0, u2=0)
+ on(3, u1=0, u2=0)
+ on(4, u1=0, u2=0);
// Solve temperature equation
solve Temperature ([dT], q)
= int2d(Th)(
Advect(du1, du2, T) * q
+ Advect(u1, u2, dT) * q
+ Graddotgrad(dT, q)
)
- int2d(Th)(
Advect(u1, u2, T) * q
+ Graddotgrad(T, q)
)
+ on(3, T = 1) // Hot bottom wall
+ on(4, T = 0); // Cold top wall
// Update fields
u1[] -= du1[];
u2[] -= du2[];
T[] -= dT[];
// Compute error with safeguard
real Lu1 = max(1e-10, u1[].linfty);
real Lu2 = max(1e-10, u2[].linfty);
real LT = max(1e-10, T[].linfty);
err = du1[].linfty / Lu1 + du2[].linfty / Lu2 + dT[].linfty / LT;
// Stop if convergence is reached
if (err < eps) break;
}
// Check convergence
if (err < eps) {
cout << "Converged after " << n << " Newton iterations.\n";
break;
}
}
```
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [February 13, 2025, 4:23pm UTC](https://community.freefem.org/t/solving-non-linear-system-of-equations-with-newton-raphson-method/3733/4 "2025-02-13T16:23:38Z")
There are several mistakes in your code:
– default labels bottom=1, right=2, top=3, left=4
– in a variational formulation you have to separate (ie in different int2d()) the terms linear with respect to the unknown and the terms which do not depend on the unknown
– boundary conditions on() : only the unknowns are allowed (du1,du2) in “velocity”, dT in “Temperature”
– it is better to resolve both velocity and temperature together in order to have fast convergence of the Newton iteration
– You have to set the boundary conditions on the updated values u1-du1, u2-du2, T-dT. Thus to set these values to 0,0,1 respectively you need to set
```auto
+ on(1,2,3,4, du1=u1, du2=u2) //(meaning u1-du1=0, u2-du2=0)
+ on(1, dT = T-1) //(meaning T-dT=1)
```
[boussinesq.edp](https://community.freefem.org/uploads/short-url/23nnM7Q1ZCVU0wp2FNfbbEjFYX5.edp) (2.6 KB)
I am doubtful of the penalty method for div u=0. Probably it gives wrong results because of no inf-sup stability. You have to use rather a variational formulation in velocity-pressure (and temperature in your case).
---
### Author:  [@Try\_everything](https://community.freefem.org/u/Try_everything)
#### Post date: [February 15, 2025, 10:47am UTC](https://community.freefem.org/t/solving-non-linear-system-of-equations-with-newton-raphson-method/3733/5 "2025-02-15T10:47:44Z")
Very thanks sir! Yes, the results are wrong. I am not aware of the stability condition. Without penalty method, I am not getting convergence. What should I look into?
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [February 15, 2025, 3:10pm UTC](https://community.freefem.org/t/solving-non-linear-system-of-equations-with-newton-raphson-method/3733/6 "2025-02-15T15:10:27Z")
The inf-sup condition of Babuska-Brezzi is very classical, see for example
> **[Ladyzhenskaya–Babuška–Brezzi condition](https://en.wikipedia.org/wiki/Ladyzhenskaya%E2%80%93Babu%C5%A1ka%E2%80%93Brezzi_condition)**
>
> In numerical partial differential equations, the Ladyzhenskaya–Babuška–Brezzi (LBB) condition is a sufficient condition for a saddle point problem to have a unique solution that depends continuously on the input data. Saddle point problems arise in the discretization of Stokes flow and in the mixed finite element discretization of Poisson's equation. For positive-definite problems, like the unmixed formulation of the Poisson equation, most discretization schemes will converge to the true solution...
see also “the system of Stokes for fluids” in FF documentation p.69
You have to use a couple of spaces for velocity and pressure satisfying the inf-sup condition. You can take for example P2 for velocity and P1 for pressure (Taylor-Hood choice). This gives
[boussinesq.edp](https://community.freefem.org/uploads/short-url/mK6vFQ3zc8h6rT1jwAhhdelgiL7.edp) (3.0 KB)