# Negative values in solution

**URL:** <https://community.freefem.org/t/negative-values-in-solution/2642>\
**Category:** General Discussion\
**Created:** [August 4, 2023, 2:26pm UTC](https://community.freefem.org/t/negative-values-in-solution/2642 "2023-08-04T14:26:52Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Anita](https://avatars.discourse-cdn.com/v4/letter/a/49beb7/32.png) [@Anita](https://community.freefem.org/u/Anita)\
**Post date:** [August 4, 2023, 2:26pm UTC](https://community.freefem.org/t/negative-values-in-solution/2642/1 "2023-08-04T14:26:52Z")

</div>

I’m solving equations describing density of some physical tissues so obviously negative values have no sense in real world however I get some during calculation (and it’s getting worse with each time step).  
Have you faced such issue? What kind of clever techniques can be used not to influence resolution too much. I’ve initial and boundary condition, already decreased time step and no further clues.

---

<div class="post-metadata">

**Author:** ![jbl](https://avatars.discourse-cdn.com/v4/letter/j/e47c2d/32.png) [@jbl](https://community.freefem.org/u/jbl)\
**Post date:** [August 4, 2023, 2:40pm UTC](https://community.freefem.org/t/negative-values-in-solution/2642/2 "2023-08-04T14:40:57Z")

</div>

Hello, can you post some code?

You can check the sign of your linear/bilinear term. Also check that you’re writing your boundary condition correctly (the sign of the normal may have an influence).

---

<div class="post-metadata">

**Author:** ![Anita](https://avatars.discourse-cdn.com/v4/letter/a/49beb7/32.png) [@Anita](https://community.freefem.org/u/Anita)\
**Post date:** [August 4, 2023, 6:52pm UTC](https://community.freefem.org/t/negative-values-in-solution/2642/3 "2023-08-04T18:52:51Z")

</div>

Where should I check the sign? In problem definition?

Don’t want to place the whole code (it’s long). I have system of 3 equations, below 2 most connected. All functions, macros, constants are defined. Below problem is solved in 5-10 steps but etaa become negative.

problem Prob([etai,etaa],[o1,o3],solver=LU)=  
//1st  
int2d(Th)(etai/dt_o1) -int2d(Th)(etaiold/dt_o1)  
-int2d(Th)(Detai\*(dx(etai)\*dx(o1) + dy(etai)_dy(o1)))  
-int2d(Th)(wektoru(uv1, uv2)'grad(etai)o1)  
+int2d(Th)(phi(x,y)etaio1)  
//3rd  
+int2d(Th)(etaa/dto3) -int2d(Th)(etaaold/dto3)  
-int2d(Th)(Detaa_(dx(etaa)_dx(o3) + dy(etaa)dy(o3)))  
-int2d(Th)(wektoru(uv1, uv2)'grad(etaa)o3)  
-int2d(Th)(phi(x,y)etaio3)  
+int2d(Th)(alphaetaao3)  
+int2d(Th)(omegaaetaa_o3)

//boundary conditions  
+on(C04, etai = w2, etaa = 0)

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [August 5, 2023, 11:06am UTC](https://community.freefem.org/t/negative-values-in-solution/2642/4 "2023-08-05T11:06:03Z")

</div>

Can you upload the complete code? Probably the best thing to do is  
work through the discretezation eqns and see how the problem  
occurs with coarse time or space steps. Sign mistakes are common  
too ( esp for me doing int by parts in head lol ). An alterantive is  
to work with logs but that is not easy either.

---

<div class="post-metadata">

**Author:** ![Anita](https://avatars.discourse-cdn.com/v4/letter/a/49beb7/32.png) [@Anita](https://community.freefem.org/u/Anita)\
**Post date:** [August 6, 2023, 4:25pm UTC](https://community.freefem.org/t/negative-values-in-solution/2642/5 "2023-08-06T16:25:56Z")

</div>

Slightly simplified code below but the main parts are there. Maybe I’m making some silly mistake

border C01(t=0.06, 2.08){x=t; y=0.002;}  
border C02(t=0.002, 1.002){x=2.08; y=t;}  
border C03(t=2.08, 0.06){x=t; y=1.002;}  
border C04(t=1.002, 0.002){x=0.06; y=t;}  
int n = 100;  
mesh Th = buildmesh(C01(n) + C02(n/2) + C03(n) + C04(n/2));

fespace Vh1(Th,P2);  
// Parameters  
real Dgamma = 1000;  
real Detai = 20;  
real Detaa = 40;  
real gammahypo = 10;  
real fi0 = 0.9;

real w1=60;  
real w2=50 ;  
real alpha=0.5;  
real beta=0.85;  
real gamma0=0.85;  
real omegaa=0.03010299956;  
real t=0, T=5, dt=1;

Vh1 etai, etaa, gamma,utest,  
etaiold, etaaold, gammaold,  
o1,o2,o3,  
uv1 = 2,  
uv2 = 2;

//initial conditions  
etai = 0;  
gamma = gamma0;  
etaa = 0;

//functions  
macro grad(ff) [dx(ff),dy(ff)] //EOM  
macro wektoru(w1,w2) [w1,w2] //EOM

func real fi(real xx, real yy)  
{if (gamma(xx,yy) \< gammahypo)  
{return fi0;}  
else {return 0; }  
}

problem Prob([etai,gamma, etaa],[o1,o2,o3],solver=LU)=  
int2d(Th)(etai/dt_o1) -int2d(Th)(etaiold/dt_o1)  
-int2d(Th)(Detai\*(dx(etai)_dx(o1) + dy(etai)dy(o1)))  
-int2d(Th)(wektoru(uv1, uv2)'grad(etai)o1)  
+int2d(Th)(fi(x,y)etaio1)  
//2nd  
+int2d(Th)(gamma/dto2) -int2d(Th)(gammaold/dto2)  
-int2d(Th)(Dgamma(dx(gamma)dx(o2) + dy(gamma)dy(o2)))  
-int2d(Th)(wektoru(uv1, uv2)'grad(gamma)o2)  
+int2d(Th)(betagammao2)  
//3rd  
+int2d(Th)(etaa/dto3) -int2d(Th)(etaaold/dto3)  
-int2d(Th)(Detaa_(dx(etaa)_dx(o3) + dy(etaa)dy(o3)))  
-int2d(Th)(wektoru(uv1, uv2)'grad(etaa)o3)  
-int2d(Th)(fi(x,y)etaio3)  
+int2d(Th)(alphaetaao3)  
+int2d(Th)(omegaaetaa_o3)  
//boundary conditions  
+on(C04, etai = w2, gamma = w1, etaa = 0)  
;

for(real t=t;t\<=T;t+=dt)  
{

cout\<\<"===================================Solution in step t: "\<\<t\<\<endl;

Prob;  
//plot(etai,value=1,fill=1,wait=0, cmm="solution etai at time: "+t);  
//plot(gamma,value=1,fill=1,wait=0, cmm="solution gamma at time: "+t);  
//plot(etaa,value=1,fill=1,wait=0, cmm="solution etaa at time: "+t);

etaiold=etai; etaaold=etaa; gammaold=gamma;

}

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [August 6, 2023, 5:33pm UTC](https://community.freefem.org/t/negative-values-in-solution/2642/6 "2023-08-06T17:33:13Z")

</div>

It is not so simple to find a numerical method such that the solution is positive.with finite element method.

1. for P2 finite element no method exist to day.
2. for P1 finite element it is generally ok with mass lump quadrature if the mesh is acute but it is hard to build a mesh with only acute triangle, but in general if the no acute triangles are not in bad place it works.

for mass lumping for P1 add

change

```auto
fespace Vh1(Th,P2);
...
int2d(Th)(etai/dto1) -int2d(Th)(etaiold/dto1)

```

in

```auto
fespace Vh1(Th,P1);
....
int2d(Th,qft=qf1pTlump)(etai/dto1) -int2d(Th,qft=qf1pTlump)(etaiold/dto1)

```

---

<div class="post-metadata">

**Author:** ![Anita](https://avatars.discourse-cdn.com/v4/letter/a/49beb7/32.png) [@Anita](https://community.freefem.org/u/Anita)\
**Post date:** [August 7, 2023, 6:02pm UTC](https://community.freefem.org/t/negative-values-in-solution/2642/7 "2023-08-07T18:02:21Z")

</div>

Thanks for answering.  
I introduced suggested modifications, results changed but still getting negative values.