有限元解一维poisson方程 matlab

一维Poisson方程的一般形式为：

$$
\frac{\partial^2 u(x)}{\partial x^2} = f(x)
$$

其中，$u(x)$是未知的函数，$f(x)$是已知的函数。

使用有限元方法求解一维Poisson方程的步骤为：

1. 离散化区域

将求解区域离散化成若干个网格，对于一维情况，通常使用等距网格。假设求解区域为$[0, L]$，将其等分为$n$个网格，每个网格长度为$h=\frac{L}{n}$。

2. 建立有限元空间

在每个网格上建立有限元空间，通常使用线性元，即局部插值函数为：

$$
N_1(x)=\frac{x_{j+1}-x}{h},\quad N_2(x)=\frac{x-x_j}{h}
$$

其中，$x_j$和$x_{j+1}$分别是第$j$个网格的左右端点。

3. 组装刚度矩阵和载荷向量

根据有限元方法的基本思想，将每个网格的刚度矩阵和载荷向量组装成整个区域的刚度矩阵和载荷向量。对于一维Poisson方程，每个网格的刚度矩阵为：

$$
K_j = \frac{1}{h} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}
$$

每个网格的载荷向量为：

$$
f_j = \begin{bmatrix} \int_{x_j}^{x_{j+1}}f(x)N_1(x)dx \\ \int_{x_j}^{x_{j+1}}f(x)N_2(x)dx \end{bmatrix}
$$

4. 边界条件处理

根据边界条件，将整个区域的刚度矩阵和载荷向量进行修正。对于一维Poisson方程，通常有三种边界条件：Dirichlet边界条件、Neumann边界条件和Robin边界条件。

- Dirichlet边界条件：$u(a)=\alpha$和$u(b)=\beta$，其中$a$和$b$分别是求解区域的左右端点，$\alpha$和$\beta$是已知的常数。将刚度矩阵和载荷向量的第一行和第$n$行全部清零，然后将对角线上的两个元素分别设为1，将载荷向量中对应的元素设为边界条件的值即可。
- Neumann边界条件：$\frac{\partial u}{\partial x}(a)=\alpha$和$\frac{\partial u}{\partial x}(b)=\beta$，其中$a$和$b$分别是求解区域的左右端点，$\alpha$和$\beta$是已知的常数。将载荷向量的第一和最后一个元素分别加上$h\alpha$和$h\beta$即可。
- Robin边界条件：$u(a)+\alpha\frac{\partial u}{\partial x}(a)=\beta$和$u(b)+\gamma\frac{\partial u}{\partial x}(b)=\delta$，其中$a$和$b$分别是求解区域的左右端点，$\alpha$、$\beta$、$\gamma$和$\delta$是已知的常数。类似于Dirichlet边界条件，将刚度矩阵和载荷向量的第一行和第$n$行全部清零，然后将对角线上的两个元素分别设为$1+\alpha h$和$1+\gamma h$，将载荷向量中对应的元素设为边界条件的值即可。

5. 求解线性方程组

将修正后的刚度矩阵和载荷向量带入线性方程组$Ku=f$中，使用MATLAB中的求解器（如$\backslash$或者inv函数）求解出未知的系数向量$u$。

下面是MATLAB代码实现：

% 定义求解区域和网格数
L = 1; n = 10; h = L/n;

% 定义有限元空间
phi1 = @(x) (x(2)-x)/(h);
phi2 = @(x) (x-x(1))/(h);

% 组装刚度矩阵和载荷向量
K = zeros(n,n);
f = zeros(n,1);
for j = 1:n-1
    Kj = 1/h*[1,-1;-1,1];
    fj = [integral(@(x) f(x).*phi1(x+j*h-[h,0]),j*h-h,j*h), ...
          integral(@(x) f(x).*phi2(x+j*h-[0,h]),j*h,j*h+h)]';
    K(j:j+1,j:j+1) = K(j:j+1,j:j+1) + Kj;
    f(j:j+1) = f(j:j+1) + fj;
end

% 处理边界条件
% Dirichlet边界条件
% K(1,:) = 0; K(1,1) = 1; f(1) = u0;
% K(n,:) = 0; K(n,n) = 1; f(n) = un;

% Neumann边界条件
% f(1) = f(1) + h*alpha;
% f(n) = f(n) + h*beta;

% Robin边界条件
K(1,:) = 0; K(1,1) = 1+alpha*h; f(1) = beta;
K(n,:) = 0; K(n,n) = 1+gamma*h; f(n) = delta-gamma*h*u(n-1);

% 解线性方程组
u = K\f;

% 绘制数值解
x = linspace(0,L,n+1);
plot(x,[0;u;0],'o-');