numerical-methods - Freefem++:用数值函数求解泊松方程

Question

我正在使用 Freefem++ 来求解泊松方程

梯度^2 u(x,y,z) = -f(x,y,z)

当我有 f 的解析表达式时，它效果很好，但现在我有一个数字定义的 f(即在网格上定义的一组数据)，我想知道是否仍然可以使用 Freefem++。

即典型代码(对于本例中的二维问题)如下所示

mesh Sh= square(10,10); // mesh generation of a square
fespace Vh(Sh,P1); // space of P1 Finite Elements
Vh u,v; // u and v belongs to Vh

func f=cos(x)*y; // analytical function

problem Poisson(u,v)= // Definition of the problem
    int2d(Sh)(dx(u)*dx(v)+dy(u)*dy(v)) // bilinear form
    -int2d(Sh)(f*v) // linear form
    +on(1,2,3,4,u=0); // Dirichlet Conditions
Poisson; // Solve Poisson Equation
plot(u); // Plot the result

我想知道是否可以用数值方式而不是分析方式来定义 f。

Answer 1

网格和空间定义

我们定义一个 Nx=10 网格和 Ny=10 的正方形单元，这在 x 轴上提供了 11 个节点，y 轴上也有相同的节点。

int Nx=10,Ny=10;
int Lx=1,Ly=1;
mesh Sh= square(Nx,Ny,[Lx*x,Ly*y]); //this is the same as square(10,10)
fespace Vh(Sh,P1); // a space of P1 Finite Elements to use for u definition

条件和问题陈述

我们不会使用求解，但我们将处理矩阵(一种使用 FreeFem 求解的更复杂的方法)。

首先，我们为我们的问题定义 CL(狄利克雷问题)。

varf CL(u,psi)=on(1,2,3,4,u=0); //you can eliminate border according to your problem state
    Vh u=0;u[]=CL(0,Vh);
    matrix GD=CL(Vh,Vh);

然后我们定义问题。我建议使用宏，而不是编写 dx(u)*dx(v)+dy(u)*dy(v) ，因此我们将 div 定义如下，但要注意宏以 结束// 不是 ;。

 macro div(u) (dx(u[0])+dy(u[1])) //

因此泊松双线性形式变为:

varf Poisson(u,v)= int2d(Sh)(div(u)*div(v));

提取刚度矩阵后

matrix K=Poisson(Vh,Vh);
matrix KD=K+GD; //we add CL defined above

我们继续求解，UMFPACK是FreeFem中的求解器，对此没有太多关注。

set(KD,solver=UMFPACK);

这就是您所需要的。您想要在某些特定节点上定义函数 f 的值。我将告诉你这个 secret :泊松线性形式。

real[int] b=Poisson(0,Vh);

您可以在任何您想要执行的节点定义函数 f 的值。

b[100]+=20; //for example at node 100 we want that f equals to 20
b[50]+=50; //and at node 50 , f equals to 50 

我们解决了我们的系统。

u[]=KD^-1*b; 

终于我们得到了情节。

plot(u,wait=1);

希望这对你有帮助，感谢我的实习导师 Olivier，他总是给我特别在 FreeFem 上的技巧。我测试过，效果很好。祝你好运。