写一个一维偏微分方程的例题并用有限元法求解，用Matlab演示

我们考虑求解如下的一维波动方程：

$$\frac{\partial^2 u}{\partial t^2} - c^2 \frac{\partial^2 u}{\partial x^2} = f(x,t)$$

其中，$u(x,t)$表示波动的位移，$c$为波速，$f(x,t)$为外部激励力。我们将求解区域$[0,1]$分成$N$个小段，每个小段长度为$h=\frac{1}{N}$。定义$x_i = i h$为第$i$个节点的位置，$u_i = u(x_i)$为第$i$个节点的位移。则在每个小段上，我们可以将位移分布近似为线性函数：

$$u(x) \approx u_i + \frac{u_{i+1}-u_i}{h}(x-x_i)$$

将该近似代入波动方程中，得到：

$$\frac{\partial^2 u_i}{\partial t^2} - c^2 \frac{u_{i+1}-2u_i+u_{i-1}}{h^2} = f_i$$

注意边界条件，当$x=0$或$x=1$时，位移应该是固定的，即$u(0,t)=u_0, u(1,t)=u_N$。因此，我们可以将边界节点的位移固定，只对内部节点的位移进行求解。整理上式，得到：

$$\frac{d^2u_i}{dt^2} = c^2 \frac{u_{i+1}-2u_i+u_{i-1}}{h^2} + f_i$$

使用欧拉法进行时间离散化，得到：

$$\frac{u_i^{n+1}-2u_i^n+u_i^{n-1}}{\Delta t^2} = c^2 \frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{h^2} + f_i^n$$

其中，$u_i^n$表示第$n$个时间步长时第$i$个节点的位移，$\Delta t$为时间步长。将上式整理，得到：

$$u_i^{n+1} = 2u_i^n - u_i^{n-1} + \frac{c^2\Delta t^2}{h^2}(u_{i+1}^n-2u_i^n+u_{i-1}^n) + \Delta t^2 f_i^n$$

这就是一维波动方程的数值解法。下面是用matlab实现该方法的代码：

% 定义求解区间和节点数
L = 1; % 求解区间长度
N = 100; % 节点数，网格大小为 h = L / N

% 定义时间步长和求解时间段
dt = 0.001;
t = 0:dt:1;

% 定义初始位移分布和边界条件
u0 = zeros(1, N);
u0(1) = 0; % 左边界位移为0
u0(N) = 0; % 右边界位移为0

% 定义热传导系数和外部激励力
c = 1;
f = zeros(size(u0)); % 外部激励力为0

% 定义初始速度分布
v0 = zeros(1, N);
v0(2:N-1) = (u0(3:N)-u0(1:N-2))/(2*h); % 使用中心差分法求解速度

% 进行时间循环求解
u = u0;
v = v0;
for i = 2:length(t)
    u_new = zeros(size(u));
    for j = 2:N-1 % 内部节点进行计算
        u_new(j) = 2*u(j) - u(j-1) + c^2*dt^2/h^2 * (u(j+1) - 2*u(j) + u(j-1)) + dt^2*f(j);
    end
    v_new = zeros(size(v));
    for j = 2:N-1 % 内部节点进行计算
        v_new(j) = v(j) + c^2*dt^2/h^2 * (v(j+1) - 2*v(j) + v(j-1)) + dt^2/h*f(j);
    end
    u = u_new;
    v = v_new;
end

% 可视化求解结果
for i = 1:length(t)
    plot((1:N)*h, u);
    xlim([0 L]);
    ylim([-1 1]);
    xlabel('x');
    ylabel('u');
    drawnow;
end

在上述代码中，我们采用了欧拉法进行时间离散化，使用for循环对每个时间步长进行求解，最后使用plot函数将求解结果可视化出来。由于波动方程需要同时求解位移和速度，因此我们还需要计算初始速度分布，并在时间循环中同时求解位移和速度。