向大神求救，fortran中计算第一类和第二类的虚宗量贝塞尔函数I和K，库函数和代码都行！

关于虚宗量Bessel函数的Fortran代码可参见

http://jean-pierremoreaupagesperso-orangefr/f_besselhtml

第一类I

http://jean-pierremoreaupagesperso-orangefr/Fortran/tzbesi_f90txt

第一类J

http://jean-pierremoreaupagesperso-orangefr/Fortran/tzbesj_f90txt

第二类Y

http://jean-pierremoreaupagesperso-orangefr/Fortran/tzbesy_f90txt

第三类K

http://jean-pierremoreaupagesperso-orangefr/Fortran/tzbesk_f90txt

辅助的复函数库

http://jean-pierremoreaupagesperso-orangefr/Fortran/complex_f90txt

↓

文档介绍：贝塞尔函数的母函数（生成函数） 1 母函数（生成函数） 加法公式 利用母函数公式 二、球贝塞尔函数 解：以柱轴为 z 轴，底面为 平面，依题意定解问题表示为 用分离变数法求解定解问题，必须上下底、或柱侧面为齐次边界条件，为此令 由于在柱坐标系中 含 ，因此 z的一次函数 的调和量为零，这样定解转化为 令 ，代入泛定方程分离变数得 因 的定解条件与 无关，是轴对称定解问题，所以 常数， 。由上下底面齐次边界条件 ，该条件与 的微分方程构成本征值问题，其解为 当 、 时， 的解为 因为当 时， 、 ，所 以满足“ 有限”的解为 ，这样 由柱侧面的边界条件得 一、球坐标系波动、输运方程的分离变量 36 球贝塞尔方程 回顾：分离变量 1 波动方程时空变量的分离 自由波动方程有意义的解是振荡解，要求 满足的方程称为亥姆霍兹方程 2 输运方程时空变量的分离 无源输运方程有意义的解是收敛解，要求 为实数 满足的方程称为亥姆霍兹方程 3 亥姆霍兹方程在球面坐标系中的分离变数 亥姆霍兹方程 令 ，代入上式得 球函数方程 l阶球贝塞尔方程 l阶球贝塞尔方程 若令 ，代入上式 阶贝塞尔方程 其通解为： l阶球贝塞尔方程的通解为 1 定义 球贝塞尔函数： 球诺伊曼函数： 球汉克尔函数： 2球贝塞尔方程的解 球贝塞尔方程 有两个线性独立解为 和 ，通解为 三、球贝塞尔函数的基本性质 1 2 递推公式 ， 或 ， 或 例1： 3 当 时 当 时 4 球贝塞尔函数及其导数 、 有无穷多 个正零点。 例2 的正零点为 5 球贝塞尔函数正交关系与模方 设 是半径为 的球面上的齐次边界条件： 或 或 的第个 正根，则 称为球贝塞尔函数的模方。 第一类齐次边界条件： 第二类齐次边界条件： 第三类齐次边界条件： 6 以球贝塞尔函数 为基的广义傅里叶级数 在区间 上，以球贝塞尔函数 为基，

2011贝塞尔方程拉普拉斯方程在柱坐标系下的分离变量得出了一般的贝塞尔方程，由于贝塞尔方程的普遍性，我们还能从其它典型的数学物理定解问题来导出贝塞尔方程的一般形式． 考虑固定边界的圆膜振动，可以归结为下述定解问题 （2011）其中 为已知正数， 为已知函数．这个定解问题因宜于使用柱坐标，从而构成柱面问题．（由于是二维问题，即退化为极坐标） 设 ，对泛定方程分离变量（取 ）得（2012）（2012）再令 ，得到 (2013) （2014） 令 于是(2015)得到 （2015）边界条件为 方程（2015）称为 阶贝塞尔微分方程．这里 和 可以为任意数．2012 贝塞尔方程的解 通过数学物理方程的幂级数求解方法可以得出结论：（1）当 整数时，贝塞尔方程(2016)的通解为 （2016）其中 为任意常数， 定义为 阶第一类贝塞尔函数（简称为贝塞尔函数）．但是我们应该注意到：当 整数时，有 ，故上述解中的 与 是线性相关的，所以(2016)成为通解必须是 整数（2）当 取任意值时：定义第二类贝塞尔函数 ，这样贝塞尔方程的通解可表示为 （2017）（第二类贝塞尔函数也可写成 ，本书之所以选取 ，是因为它又称为诺依曼函数，取第一大写字母）下面我们会看到不管 是否为整数，上式均成立．(3) 当 取任意值时： 实际上由第一、二类贝塞尔函数还可以构成线性独立的第三类贝塞尔函数 ，又称为汉克尔函数． （2018）并分别将 称为第一种和第二种汉克尔函数． 最后，总结 阶贝塞尔方程的通解通常有下列三种形式： （i） 整数） （ii） 可以取任意数）（iii） 可以取任意数） 注明：第一、第二、第三类贝塞尔函数分别称为贝塞尔函数、诺依曼函数、汉克尔函数，还可以分别简称为第一、第二、第三类柱函数．

其实，这个问题本身并不困难，要点只在于怎样计算贝塞尔函数。

MATLAB提供了计算贝塞尔函数的函数，详情请参见另一个问题的回答。

 

首先，试图使用符号数学工具箱求解析解，代码如下：

syms x y

eq1=(cos(x)-1/3(cos(x))^3)+2(1000+683-22480)/ 

    (3(1000-683))-y(sin(x))^2/0000115-2sin(x) 

    sin(x-0959931)/0000115^2/(98(1000-683)/0071);

eq2=sin(0959931-x)besselk(0,(98(1000-683)/0071)^ 

    050000115sin(x))+(98(1000-683)/0071)^05y 

    besselk(1,(98(1000-683)/0071)^050000115sin(x));

[x,y]=solve(eq1,eq2)

得到的结果为：

x =

72408611253328779611417644360562

y =

-85870621389208280165235516006804e-6

这是因为无法求得解析解，因而调用了数值方法求解得到的结果（注意，在不同MATLAB版本中的处理方式可能存在差别，我使用Maple内核的2008a求解得到上述结果，但使用MuPad内核的2012b则得到复数解）。

 

由于这里求得的x不符合0-pi区间的要求，所以，考虑直接用数值方法求解方程（使用优化工具箱fsolve函数）。代码如下：

function zd

x0 = [pi/2; -1e-5];

options=optimset('Display','iter');

x = fsolve(@eqs, x0, options);

fprintf('x = %6g, y = %6g\n', x);

fprintf('Eq1(x, y) = %6g, Eq1(x, y) = %6g\n', eqs(x));

function f = eqs(X)

x = X(1);

y = X(2);

f(1)=(cos(x)-1/3(cos(x))^3)+2(1000+683-22480)/ 

    (3(1000-683))-y(sin(x))^2/0000115-2sin(x) 

    sin(x-0959931)/0000115^2/(98(1000-683)/0071);

f(2)=sin(0959931-x)besselk(0,(98(1000-683)/0071)^ 

    050000115sin(x))+(98(1000-683)/0071)^05y 

    besselk(1,(98(1000-683)/0071)^050000115sin(x));

得到的结果：

x = 0957676, y = -858706e-007

Eq1(x, y) = -515143e-014, Eq1(x, y) = -138778e-017

我觉得没有太多需要的说明的，所以就不多写了，有问题请追问。

 

另外，我注意到，前面用符号数学工具箱求出的解与使用fsolve得到的结果相差2pi，有兴趣可自行验证。

 

==================================

由于系统抽风，另一个问题的答案“正在提交中”，我把主要内容贴到这里，供参考。

 

MATLAB提供了计算贝塞尔函数的函数，具体包括：

besselj - 第一类贝塞尔函数，或简称贝塞尔函数；

bessely - 第二类贝塞尔函数，又称诺伊曼函数（Neumann function）；

besseli - 第一类修正贝塞尔函数；

besselk - 第二类修正贝塞尔函数；

besselh - 第三类贝塞尔函数，又称汉克尔函数（Hankel function）。

 

这几个函数的调用语法基本相同，例如

J = besselj(nu,Z)

J = besselj(nu,Z,1)

[J,ierr] = besselj(nu,Z)

其中，nu为贝塞尔函数的阶数，Z为函数自变量。阶数必须为实数，但Z可以是复数。

 

就你的问题而言，非常简单，K0(x)、K1(x)在MATLAB中的表达式分别为besselk(0,x)、besselk(1,x)。

另外值得一提的是，上述函数是MATLAB基本模块提供的特殊函数（也就是说不需要任何附加工具箱），采用数值方法计算；而符号数学工具箱则提供了第一和第二类的4个贝塞尔函数，名称和调用方式都与基本模块的4个函数完全一致，但支持微分、积分等符号运算。

由于x=0时，K1为无穷大（inf），y也为无穷大，在使用trapz函数时就认为b是无效值（nan）。解决的方法是人为把x的0值改为eps值（eps系统默认为无穷小量）。即

x=(eps:00001:1)。运行结果值，b=042596

[X,Y] = meshgrid(-4:0025:2,-15:0025:15)

H = besselh(0,1,X+iY);

contour(X,Y,abs(H),0:02:32),hold on

contour(X,Y,(180/pi)angle(H),-180:10:180);hold off

参考了《数学物理方程的MATLAB解法与可视化》（彭芳麟），p26。