Dear all
Can we write third or fourth-order derivatives in FreeFEM?
However, it is showing identifiers errors while writing third-order derivatives ( like dxxx( u ) ) .
Kindly suggest me in solving this issue.
1 Like
Applying the Greens theorem multiple times will reduce the degree of PDEs term.
It solved my problem.
What about 6th order?
You may use macro, such as:
macro grad(u) [dx(u),dy(u)]//
macro div(u1,u2) (dx(u1)+dy(u2))//
macro Vgradu(v1,v2,u) (grad(u)'*[v1,v2])//
macro Grad(u1,u2) [grad(u1),grad(u2)]//
macro VgradU(v1,v2,u1,u2) [Vgradu(v1,v2,u1),Vgradu(v1,v2,u2)]//
In 1d case suppose we want to write u_xxxxxx means in weak formulation we can pass the the derivative and we can get (u_xxx, v_xxx) or more we can pass.
Now how we can tackle this?
You can write your equation in a system form.
For example if you want to solve
u_{xxxxxx}=f,
you can write
u_x=u_1,\quad (u_1)_x=u_2,\quad (u_2)_x=u_3,\quad (u_3)_x=u_4,\quad
(u_4)_x=u_5,\quad (u_5)_x=f.
Then you have the unknowns u,u_1,u_2,u_3,u_4,u_5 in P1
and you write the variational formulation for v,v_1,v_2,v_3,v_4,v_5 in P1
\int \partial_x(u_i) v_i=\int u_{i+1}v_i,\quad i=0,...,4,
\int \partial_x(u_5)v_5=\int fv_5.
In order to have a well-posed system you need to introduce boundary conditions, by adding boundary terms.
A 2d example is to solve the Bilaplace problem \Delta^2 u=f=1 in the unit disc with boundary conditions u=0 and \nabla u=0 on the boundary.
bilaplacian.edp (1.0 KB)
In this file two methods are used. The first is by the Morley finite element space. The second is by the above method, where you have the unknowns
ub,
uxb=dx(ub)
uyb=dy(ub)
ulb=dx(uxb)+dy(uyb)
Dear sir ,
sorry for the inconvenience, i couldn’t understand it properly.
can you please explain it with 1d case ?
Thank you.
As an example we can consider the 1d problem \partial^4_x u=1 in (0,1) with boundary condition u=0 and \partial_x u=0 on the whole boundary, which is the bilaplacian problem in 1d. The solution is given by u=x^2(1-x)^2/24.
You can solve it by introducing the unknowns
u, u1=dx(u), u2=dx(u1)
It gives
solve bilap1d(u,u1,u2,v,v1,v2)=
-int1d(Th)(dx(u2)*dx(v2))
+int0d(Th)(dx(u2)*v2*N.x)
-int1d(Th)(v2)
+int1d(Th)((dx(u)-u1)*v)
+int1d(Th)(-u1*dx(v1)-u2*v1)
+int0d(Th)(1.e30*u*v2)
;
Whole file is
bilap1d.edp (818 Bytes)
It contains also another possible resolution with only two unknowns u,u2
1 Like
Dear sir,
I am very new to freefem . So in the context, bilap1d(u,u1,u2,v,v1,v2) , if we say u1 means does it consider it as dx(u) and so on? and can you please explain about the code you provided ?
Thank you.
The variational formulation that is solved by bilap1d() is: find u, u_1, u_2 in P1 such that for all v, v_1, v_2 in P1
-\int \partial_x u_2\partial_x v_2+\int_\Gamma \partial_x u_2 v_2 N-\int v_2
+\int (v\partial_x u - v u_1 )
+\int(-u_1\partial_x v_1 - u_2 v_1)+\int_\Gamma \xi u v_2\ =0,
where \int stands for \int_0^1,
\int_\Gamma f=f(0)+f(1), and N is the exterior normal i.e. N(0)=-1, N(1)=1,
and \xi=10^{30} is a large coefficient.
In this formulation, you can take successively v_1=0, v_2=0, v arbitrary, then
v=0, v_2=0, v_1 arbitrary, and finally v=0, v_1=0, v_2 arbitrary. The result is that it is equivalent to write
-\int \partial_x u_2\partial_x v_2+\int_\Gamma \partial_x u_2 v_2 N-\int v_2+\int_\Gamma \xi u v_2\ =0 for any v_2 in P1,
\int (v\partial_x u - v u_1 )=0 for any v in P1,
\int(-u_1\partial_x v_1 - u_2 v_1)=0 for any v_1 in P1.
The second line is an approximation of \partial_x u=u_1.
The third line is an approximation of \partial_x u_1=u_2 and u_1(0)=0, u_1(1)=0, because
\int -u_1\partial_x v_1=\int v_1\partial_x u_1-\int_\Gamma u_1 v_1.
The first line is an approximation of \partial^2_{xx}u_2=1 and u(0)=0, u(1)=0, because
-\int \partial_x u_2\partial_x v_2+\int_\Gamma \partial_x u_2 v_2 N=\int v_2\partial^2_{xx}u_2 (but notice that \partial^2_{xx} u_2 is not a function, it contains Dirac masses).
Another possibility which is closer to the formulation of my first message (and maybe easier to understand) is to use the unknowns u, u1=dx(u), u2=dx(u1), u3=dx(u2). This gives
fespace V(Th,P1);
V u,u1,u2,u3,v,v1,v2,v3;
solve bilap1d(u,u1,u2,u3,v,v1,v2,v3)=
int1d(Th)((dx(u2)-u3)*v2)
+int1d(Th)(dx(u3)*v3)
-int1d(Th)(v3)
+int1d(Th)((dx(u)-u1)*v)
+int1d(Th)(-u1*dx(v1)-u2*v1)
+int0d(Th)(1.e30*u*v3)
;
Sir , Thank you so much .
one doubt sir, for the given equation \partial^4_xu=1 how you are getting that weak formulation? and \partial_xu2=\partial^3_xu?
Thank you.
At the continuous level you just replace \partial^4_x u=1 by the set of equations
\partial_x u= u_1, \ \partial_x u_1=u_2, \ \partial_x u_2=u_3, \ \partial_x u_3=1.
The first equation gives u_1=\partial_x u, thus replacing u_1 by this value in the second equation gives u_2=\partial^2_x u. Then replacing u_2 by this value in the third equation gives u_3=\partial^3_x u. Finally replacing u_3 by this value in the fourth equation gives \partial^4_x u=1.
Thank you so much sir.
Suppose if we are having u_t,u_x and u may be some non linear term f(u) also means how to deal with that sir ?
Thank you.
For time dependence and nonlinearities in principle you can use standard approaches.
In particular for an evolution problem with u_t and a Cauchy initial data you need a time scheme. The most simple is to use implicit Euler. For nonlinearities you have to solve the problem by an iterative algorithm like the Newton method. The issue however is stability, and that depends strongly on the equation you want to solve and on the scheme you use.