计算科学 - 贝塞尔 EVP 和 fzero - 吾爱随笔录

贝塞尔 EVP 和 fzero

计算科学 matlab 非线性方程特殊功能根

2021-12-16 01:57:44

我正在尝试解决特征值问题微分方程是贝塞尔方程。解由第一类贝塞尔函数给出。第 2 类贝塞尔函数因其在处的奇异性而被省略。特征值必须使用处的边界条件来确定。使用恒等式处的边界条件可以改写为

x^{2} y^{″} + x y^{'} + x^{2} y = λ^{2} y, x \in (0, 1), y (0) = 0, y^{'} (1) = y (1)

$x^2 y''+ x y' + x^2 y = \lambda^2 y\,,\quad x\in(0,1)\,,\quad y(0)=0\,,\quad y'(1)=y(1)$

y (x) = J_{λ} (x)

$y(x)=J_\lambda(x)$

x = 0

$x=0$

λ

$\lambda$

x = 1

$x=1$

J_{λ}^{'} (x) = \frac{1}{2} (J_{λ - 1} (x) - J_{λ + 1} (x))

$J'_\lambda(x)=\frac{1}{2}\left(J_{\lambda-1}\left(x\right)-J_{\lambda+1}\left(x\right)\right)$

x = 1

$x=1$

J_{λ - 1} (1) - 2 J_{λ} (1) - J_{λ + 1} (1) = 0 .

$J_{\lambda-1}\left(1\right)-2J_{\lambda}\left(1\right)-J_{\lambda+1}\left(1\right)=0\,.$ 现在我想使用 MATLAB/GNU Octave 的 fzero 例程用上面的等式确定

λ

$\lambda$ 代码如下。

%% initial guess
xguess=1;
%% function handle
f=@(l) besselj(l-1,1)-2*besselj(l,1)+besselj(l+1,1)
%% set tolerance
opts=optimset('TolX',1e-12);
%% determine zero
xzero=fzero(f,xguess,opts);

如果您绘制函数 f, click here，您会发现大多数零都在 $-\infty<x<0$ 中。

要获得排序序列 $\lambda_{k+1}<\lambda_k$ ，我想仅在负 $x$ 方向上从“最后一个”零搜索到“下一个”零。有没有办法在 fzero 例程中实现这一点？

预先感谢！

2个回答

正如您在, $\log|f(\lambda)|$

f (λ) = J_{λ - 1} (1) - 2 J_{λ} (1) - J_{λ + 1} (1),

$f(\lambda) = J_{\lambda-1}(1) - 2J_\lambda(1) -J_{\lambda+1}(1),$

在此处输入图像描述

根是非常规则的，并且大约等于（从开始，是一个例外）。 $\lambda_k$ $-k$ $k\geq0$ $\lambda_0=1.23219$

因此，获得第个特征值（）的方法是将根括号并在该区间内使用求解器（它接受括号区间作为初始猜测并保持在该区间内）。 $k$ $k\geq1$ $[-k-\frac12,-k+\frac12]$ fzero

请注意，函数的大小增长非常快（类似于），因此一旦函数变得非常大（对于大），使用可能会更好作为一个近似值。 $(2\lambda/e)^\lambda\sin\pi\lambda$ $k$ $\lambda_k=-k$

这是我现在使用的代码。

kmax=5;
lambda=zeros(kmax,1);
%% function handle
f=@(l) besselj(l-1,1)-2*besselj(l,1)-besselj(l+1,1);
%% set tolerance
opts=optimset('TolX',1e-12);
%% compute first positive eigenvalue
lambda(1)=fzero(f,[0.5,1.5],opts);
%% compute negative eigenvalues
for k=1:(kmax-1)
  [lambda(k+1),~,fl]=fzero(f,[-k-0.5,-k+0.5],opts);
  if not(fl==1)
    error('iteration    did     not     convert.')
  end
end

顺便说一句，我想用它来验证我的有限差分法。

其它你可能感兴趣的问题

上一篇导出向量化操作转置的简单方法下一篇如何设置激波管问题，以便解决方案包括具有指定马赫数的激波