在matlab中求解一维波动方程的特征问题

Question

我构建了以下代码来求解一维波动方程作为半径 r 从 0 到 pi 的函数。该等式如下图所示：

我已经将空间导数 c 前面的常数取为 1，但我通常对此进行了编码，因为我希望最终使它成为一个依赖于 r 的变量。通过假设平面波解，解可以简化为：

我现在使用有限差分并将右侧的二阶 DE 离散为矩阵 A，以尝试将问题简化为求解特征值和特征向量，即：

我生成的代码如下所示：

%eigensolver code

clear
clc

%NOTE: place % in front of Neumann code when solving using Dirchlet
%condition and vise versa so that the code only inputs one type of boundary
%condition

%Specify initial variables
ri=0;               %inital radius             
rf=pi;               %final radius
Di=0;               %initial Direchlet point
Df=0;               %final Direchlet point
Ni=0;               %initial Neumann point
Nf=0;               %final Neumann point




Nr=5;
h=abs(rf-ri)/(Nr-1);  %step size
r=ri:h:rf;
A=zeros(Nr,Nr);
b=ones(size(r));
c=b;
%code inital/final point assuming Dirichlet condition
%initial
A(1,1)=-(2*(c(1)^2))/(h^2);
A(1,2)=(c(1)^2)/(h^2);
%final
A(Nr,Nr)=-(2*(c(Nr)^2))/(h^2);
A(Nr,Nr-1)=(c(Nr)^2)/(h^2);

%code initial and final points assuming Neumann condition
%initial
%A(1,1)=-(3.*(c(1)^2))./(2*h);
%A(1,2)=(2.*(c(1)^2))./h;
%A(1,3)=-(c(1)^2)/(2*h);
%final
%A(Nr,Nr-2)=(c(end)^2)/(2*h);
%A(Nr,Nr-1)=-(2*(c(end)^2))/h;
%A(Nr,Nr)=(3*(c(end)^2))/(2*h);

%interior points
for k=2:Nr-1
    A(k,k-1)=(c(k)^2)/(h^2);
    A(k,k)=-(2*(c(k)^2))/(h^2);
    A(k,k+1)=(c(k)^2)/(h^2);
end

[evec, eigval]=eig(A)

我的问题由几个部分组成：

1) 我如何知道我是否找到了数字问题的正确答案？我可以调用matlab中的一些预先指定的函数吗？我对特征问题不是很精通，所以我不知道如何解释我的答案。

2) 如果您查看上面的代码，我可以将 Dirichlet 边界条件更改为 Neumann 边界条件。这部分代码是为这个一维案例实现这些条件的正确方法吗？

Answer 1

你可以用 PDE 来思考一个特征问题，如下所示。

如果A是一个有限维矩阵，v是它的特征向量之一，l是对应的特征值，那么：

Av = lv

现在，您可以使用 PDE 的（无限维）微分算子，而不是有限维矩阵 A。在这种情况下，您的特征向量现在是特征函数，因为 D 在函数空间上工作（希望这有助于理解这一点）。

关于您的代码：您不应在刚度矩阵中包含边界条件。你想解决离散化的问题

斧头=b，

其中 A 是您的刚度矩阵，x 是您的网格点，b 是您的网格点处的值。因此，您应该在 b 内实现边界条件。

如果你想进一步阅读，解释：

以下链接是 TUM 的数字 PDE 1 课程 Numerical methods for PDE。 在这里你会找到一些关于finite differences method的解释。 课程网页将为您提供讲义、代码示例等。