# Getting issues in a simple problem that is based on both elliptic and parabolic type equation **URL:** \ **Category:** General Discussion\ **Created:** [March 31, 2025, 2:15am UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828 "2025-03-31T02:15:19Z")\ **Posts on this page:** 18\ **Page:** 1
**Author:** ![Monirul25](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/monirul25/32/3287_2.png) [@Monirul25](https://community.freefem.org/u/Monirul25)\ **Post date:** [March 31, 2025, 2:15am UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/1 "2025-03-31T02:15:19Z")
Dear all experts and respected professors, i am wondering why the following simple equation does not giving good results\> ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/3/3dc729195604f595ced5c4128ca303c6304f1fa1.png) Deatials pdf of Question\> [help.pdf](https://community.freefem.org/uploads/short-url/w4h6ZimE8OdqbaKnYcJZylarEQ5.pdf) (156.4 KB) ls pdf\> Here, is my codes: [Heat\_Elliptic\_combined.edp](https://community.freefem.org/uploads/short-url/c8c2PsoVxjLkGLjx4sozrnPvr56.edp) (2.7 KB) Thanks in advance!. ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [March 31, 2025, 11:10am UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/2 "2025-03-31T11:10:59Z")
There are several mistakes in your code: - You miss `t+=dt` - In order to solve the Neumann problem for z you need a right-hand side with really vanishing integral at the discrete level, and to add `eps*z*phi` see the FF documentation [https://doc.freefem.org/pdf/FreeFEM-documentation.pdf](https://doc.freefem.org/pdf/FreeFEM-documentation.pdf), note at the bottom p.218 - To have vanishing integral in the right-hand side of the Laplace equation on `z`, one possibility is to decouple the problem (at each timestep): first solve the equation on `u`, then once `u` and its average can be computed, solve the equation for `z`. However this works only if there is no dependence in `z` in the equation on `u` (this is true for your example). - If you have a full coupling (dependency in `z` in the equation on `u`) and you want to solve the problem in a coupled manner (solving `u` and `z` simultaneously), you can use the Lagrange multiplier method. This means that you take the average of `u` as additional scalar unknown. Here is your code with this Lagrange multiplier method. [Heat\_Elliptic\_combined.edp](https://community.freefem.org/uploads/short-url/65EBQJYmb3iRjQa89qcgxjlvggS.edp) (3.0 KB) ---
**Author:** ![Monirul25](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/monirul25/32/3287_2.png) [@Monirul25](https://community.freefem.org/u/Monirul25)\ **Post date:** [March 31, 2025, 11:36am UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/3 "2025-03-31T11:36:40Z")
Thanks you so much for your responses. It helps me a lot. **Can this Lagrange multiplier method is applicable if the equation is looks like the following:** ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/6/6a4c819f9aa43b7f9a1ae671347209311501fee0.png) (**Does (2.1), (2.2), (2.3) can be solve using this method sir simultaneously??)**. Many many Thanks sir!. ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [March 31, 2025, 11:51am UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/4 "2025-03-31T11:51:00Z")
A priori yes, but the stability of the scheme is quite much unclear. You have to use \bar\phi instead of \bar\phi\_0 to be sure that the right-hand side of (2.1c) has vanishing integral at the discrete level and for all times. There is a condition missing in your system to fix a possible additive constant in \xi when solving (2.1c) with Neumann BC. But it may be unnecessary since adding a constant to \mu via (2.1b) will not modify (2.1a) nor (2.1d). The most simple is to set \int\xi=0 (this is what is taken in the previous code). ---
**Author:** ![Monirul25](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/monirul25/32/3287_2.png) [@Monirul25](https://community.freefem.org/u/Monirul25)\ **Post date:** [March 31, 2025, 12:06pm UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/5 "2025-03-31T12:06:50Z")
So, i will take average of phi instead of average of phi\_{0} in the code sir??. Actually, i have to solve equation (2.1a) -(2.2c) by coupling. ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [March 31, 2025, 12:10pm UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/6 "2025-03-31T12:10:12Z")
Yes. For the continuous problem one has \int\phi(t)=\int\phi\_0, but for the numerical approximation it is not true exactly. You have to take \int\phi(t). ---
**Author:** ![Monirul25](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/monirul25/32/3287_2.png) [@Monirul25](https://community.freefem.org/u/Monirul25)\ **Post date:** [March 31, 2025, 12:10pm UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/7 "2025-03-31T12:10:59Z")
Okay thank you sir. One more thing sir that this method can work in DG setting??. (That is, i have to use method in DG setting by replacing all function space in P1dc or P2dc). ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [March 31, 2025, 12:12pm UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/8 "2025-03-31T12:12:00Z")
I don’t know. You can try. ---
**Author:** ![Monirul25](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/monirul25/32/3287_2.png) [@Monirul25](https://community.freefem.org/u/Monirul25)\ **Post date:** [March 31, 2025, 12:12pm UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/9 "2025-03-31T12:12:59Z")
Yes sir. Thank you so much❤. ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [March 31, 2025, 12:33pm UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/10 "2025-03-31T12:33:49Z")
It seems that you can somehow simplify your system by setting \mu\_1=f(\phi)-\varepsilon\Delta\phi. Then you get \phi\_t+\nabla\cdot(\phi u)-\varepsilon\Delta \mu\_1+\varepsilon\theta(\phi-\bar\phi)=0, \mu\_1=f(\phi)-\varepsilon\Delta\phi, -\Delta\xi=\theta(\phi-\bar\phi), with BC \partial\_n\phi=\partial\_n\mu\_1=\partial\_n\xi=0. In this form \xi can be computed afterwards. The coupling is only in the unknowns \phi and \mu\_1. This approach is possible as long as you solve the u separately. If you need to compute all variables \phi, \mu, \xi, u simultaneously (for stability reasons?), then the previous decoupling is not possible because \nabla \mu (and therefore \nabla\xi) is needed in (2.1d). I think that the resolution of u separately can be stable because it is involved in (2.1a) only in the advection term, which is low order with respect to the \Delta\mu term. ---
**Author:** ![Monirul25](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/monirul25/32/3287_2.png) [@Monirul25](https://community.freefem.org/u/Monirul25)\ **Post date:** [March 31, 2025, 12:46pm UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/11 "2025-03-31T12:46:22Z")
Yes sir. But one paper is there where they solved all the variables simultaneously by using pressure correction sir. They provided the stablity also. Here, is the paper\> [DG PPC of CHDS.pdf](https://community.freefem.org/uploads/short-url/rK8guUQjMAedb0epdpx38IGxO4p.pdf) (2.0 MB) ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [March 31, 2025, 12:49pm UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/12 "2025-03-31T12:49:47Z")
Ok, do what you think is best! ---
**Author:** ![Monirul25](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/monirul25/32/3287_2.png) [@Monirul25](https://community.freefem.org/u/Monirul25)\ **Post date:** [March 31, 2025, 1:08pm UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/13 "2025-03-31T13:08:38Z")
Thanks sir. DG setting is also working sir. Many many thanks for your support. ---
**Author:** ![Monirul25](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/monirul25/32/3287_2.png) [@Monirul25](https://community.freefem.org/u/Monirul25)\ **Post date:** [April 1, 2025, 7:40am UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/14 "2025-04-01T07:40:12Z")
Respected sir, i have implemented your idea of Lagrange multiplier to solve first three equations (2.1a), (2.1b), and (2.1c) in the following equations: ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/6/6a4c819f9aa43b7f9a1ae671347209311501fee0.png) Here, is a piece of code sir\>\> > Blockquote > real theta=1.; > real a=0.1; > real ac=1./a; > // Finite Element Spaces > fespace Vh(Th,P2); Vh u,w,z,phi,psi,rho; fespace Vh3(Th,[P2,P2,P2]); Vh3 [wu,ww,wz]; real uAverage; varf averagephi(u,w,z,phi,psi,rho)=int2d(Th)(phi/Th.area); real[int] avphiform=averagephi(0,Vh3); varf averagepsi(u,w,z,phi,psi,rho)=int2d(Th)(psi/Th.area); real[int] avpsiform=averagepsi(0,Vh3); varf averagerho(u,w,z,phi,psi,rho)=int2d(Th)(rho/Th.area); real[int] avrhoform=averagerho(0,Vh3); // To solve First three equations (2.1a), (2.1b), and (2.1c) simultaneously ``` varf pbuwz(u,w,z,phi,psi,rho)=int2d(Th)(u*phi/dt+a*Grad(w)'*Grad(phi)) // Due to (2.1a) -int2d(Th)(w*psi)+int2d(Th)((ac*(ukold^3-ukold)-stab*ukold)*psi)// Non-linear term for (2.1b) +int2d(Th)(stab*u*psi)//stabilization +int2d(Th)(a*Grad(u)'*Grad(psi)) // for (2.1b) +int2d(Th)(z*psi) // for equaution of z in (2.1c) -int2d(Th)(u0*v10*dx(phi)+u0*v20*dy(phi)) // advection term in (2.1a) -int2d(Th)(u0*phi/dt+f(t+dt/2.)*phi)// RHS of first equation in (2.1a) +int2d(Th)(Grad(z)'*Grad(rho)) // for equation of z in (2.1c) +int2d(Th)(1.e-8*z*rho) // Due to (2.1c) -int2d(Th)(u*rho); // Due to (2.1c) // solve matrix A=pbuwz(Vh3,Vh3,Vh3); real[int] rhs=pbuwz(0,0,Vh3); rhs=-rhs; matrix B=[[A,avrhoform],[avphiform',-1.],[avpsiform',-1.]]; set(B,solver="UMFPACK"); real[int] rhsB(rhs.n+1); rhsB=[rhs,0., 0.]; real[int] solB=B^-1*rhsB; [wu[],ww[],uAverage]=solB; u=wu; z=wz; w=ww; ``` > Blockquote > Here, u0, [v10, v20](initial velocity) are knows due to initial datas. I have used the continuous setting of the paper\>\> [DG PPC of CHDS.pdf](https://community.freefem.org/uploads/short-url/rK8guUQjMAedb0epdpx38IGxO4p.pdf) (2.0 MB) I have not given nonlinear loop and time loop here sir (ukold comes for nonlinearity) (as i know you know better than me sir). I have tried to solve as in the paper in continuous setting: ![Help](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/2/2243e7f0c8f2c4a81fd4b91c8487c99d401a0405.jpeg) But some errors are coming in matrix A. Here, is the errors\>\> ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/4/42e29842b4387c99597826446c509aa210a0141e.png) I will be very grateful sir if i get some help on this. Thanks in advance sir!. ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [April 1, 2025, 10:11am UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/15 "2025-04-01T10:11:35Z")
It has to be ```auto matrix A=pbuwz(Vh3,Vh3); real[int] rhs=pbuwz(0,Vh3); rhs=-rhs; matrix B=[[A,avrhoform],[avphiform',-1.]]; set(B,solver="UMFPACK"); real[int] rhsB(rhs.n+1); rhsB=[rhs,0.]; real[int] solB=B^-1*rhsB; [wu[],uAverage]=solB; u=wu; w=ww; z=wz; ``` ---
**Author:** ![Monirul25](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/monirul25/32/3287_2.png) [@Monirul25](https://community.freefem.org/u/Monirul25)\ **Post date:** [April 1, 2025, 10:45am UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/16 "2025-04-01T10:45:38Z")
Thank you so much sir for your great help. Now, code is working fine sir. \*\*But Can you explain in short why \*\* **matrix A=pbuwz(Vh3,Vh3) instead of matrix A=pbuwz(Vh3,Vh3,Vh3);** **still there were three unkowns in the bilinear form of pbuwz. I am not getting this sir**??. ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [April 1, 2025, 11:23am UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/17 "2025-04-01T11:23:29Z")
The declaration is `varf pbuwz(u,w,z,phi,psi,rho)` Then `A=pbuwz(Vh3,Vh3)` means that: the first Vh3 is for the set of unknowns (u,w,z)\in Vh3, the second Vh3 is for the test functions (phi,psi,rho)\in Vh3 ---
**Author:** ![Monirul25](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/monirul25/32/3287_2.png) [@Monirul25](https://community.freefem.org/u/Monirul25)\ **Post date:** [April 1, 2025, 12:00pm UTC](https://community.freefem.org/t/getting-issues-in-a-simple-problem-that-is-based-on-both-elliptic-and-parabolic-type-equation/3828/18 "2025-04-01T12:00:36Z")
Okay sir. Very nice explanation.