计算科学 - 波动方程，波从纽曼边界反弹 - 吾爱随笔录

波动方程。混合公元前。

应用 Neuman 边界条件 ( $\frac{\partial u}{\partial x}\big|_N=0$ ) 到域的 RHS

您可能会在所附图像中观察到字符串中间的锋利边缘。我确实理解它在 Dirichlet BC 反弹，但由于某种原因，它也从右侧的 Neuman 边界反弹。所以它从 L 到 R 来回反弹。

PS也试过 $u_N=u_{N-1}$ ，同样的结果。（我的节点索引是 1 到 101，即最后一个节点是101=n+1）。

您能否链接任何克服此问题的资源？或者解释一下我在边界做错了什么。

MatLab 的代码附在下面。

clc;clear all;clf; close all;
L = 1;
n = 100;
dx = L / n;
ld = 30; 
dt = sqrt(dx^2/ld);
x = 0:dx:L;
u_prev = sin(2*pi*x)+5*sin(3*pi*x);
u = dt*(3*sin(5*pi*x)) + u_prev;
u_new = u + 1;
vareps = 10^-8;
delta = 1; 
ctr = 0;
mplot = 1;
while delta > vareps
    for i=2:n
        u_new(i) = ld*dt*dt*(u(i-1)-2*u(i)+u(i+1))/dx/dx + 2*u(i) - u_prev(i);
    end
    u_new(1)=0;
    u_new(n+1) = -(-4*u_new(n)+u_new(n-1))/3; % Neuman condition du/dx = 0 at Right boundary
    if mplot > 0
        plot(x,u_new); 
        axis([0 1 -50 50])
        drawnow;
        mplot=0;
    end
    delta = max(abs(u-u_new));
    u_prev = u;
    u = u_new;
    
    mplot = mplot + delta;
    ctr = ctr + 1;
end

更新。尝试了 Thomas 方法，结果相同... :(

clc;clear;clf; close all; format long;
n = 101;
dx = 1/(n-1);
dt = 10^-4;
alpha = 3;
u0 = 0;
ux = 0;
m = 1;
yps = 10^-5;
x = 0:dx:1;
prev = ut0(x);
curr = upt0(x)*dt+prev;
next = prev + 1;
prev(1)=0; prev(n)=prev(n-1);
curr(1)=0; curr(n)=curr(n-1);
next(1)=0; next(n)=next(n-1);
a = curr;
b = curr;
ctr = 0;
A = (-1*alpha^2)/(dx^2);
B = (2*alpha^2/dx^2)+(1/dt^2);
C = -1*alpha^2/(dx^2);
while m > yps
    a(2) = 0; b(2) = u0;
    D = (2*curr-prev)/(dt^2);
    %forw
    for i = 2:n-1
        a(i+1) = (-A)/(B + C*a(i));
        b(i+1) = (D(i)-C*b(i))/(B+C*a(i));
    end
    %backw
    next(n)=b(n)/(1-a(n)); % Neuman from P_{n-1}=a_n P_n + b_n P_{n-1}=P_n
    for i = n-1:-1:1  
        next(i)=next(i+1)*a(i+1)+b(i+1);
    end
    if mod(ctr,50)==0 
        plot(x,next,'b');
        axis([0 1 -8 8]);
        drawnow;
    end
    prev  = curr;
    curr = next;
    m = max(abs(prev-curr));
    next = curr+1;
    ctr = ctr+1;
end

function y = ut0(x)
    y = sin(2*pi*x)+5*sin(3*pi*x);
end

function y = upt0(x)
    y = 3*sin(5*pi*x);
end