2个关于xy,并含修正贝塞尔函数的方程,怎样用matlab求解xy,

2个关于xy,并含修正贝塞尔函数的方程,怎样用matlab求解xy,
两个含xy的方程分别为
(cos(x)-1./3.*(cos(x)).^3)+2.*(1000+683-2.*2480)./(3.*(1000-683))-y.*(sin(x)).^2./0.000115-2.*sin(x).*sin(x-0.959931)./0.000115.^2./(9.8.*(1000-683)./0.071)=0
sin(0.959931-x).*K0((9.8.*(1000-683)./0.071).^0.5.*0.000115.*sin(x))+(9.8.*(1000-683)./0.071).^0.5.*y.*K1((9.8.*(1000-683)./0.071).^0.5.*0.000115.*sin(x))=0
x区间 0,pi
y区间 -0.000100,0
第二个方程含有修正贝塞尔函数
其中K0表示第二类零阶修正贝塞尔函数
K1表示第二类一阶修正贝塞尔函数
这里还有一个刚提问的相关问题,关于第二类零阶修正贝塞尔函数的
网址贴不上去,question/574166458.html
到时积分一并奉上

其实,这个问题本身并不困难,要点只在于怎样计算贝塞尔函数.
MATLAB提供了计算贝塞尔函数的函数,详情请参见另一个问题的回答.
 
首先,试图使用符号数学工具箱求解析解,代码如下：
syms x y
eq1=(cos(x)-1/3*(cos(x))^3)+2*(1000+683-2*2480)/ ...
    (3*(1000-683))-y*(sin(x))^2/0.000115-2*sin(x)* ...
    sin(x-0.959931)/0.000115^2/(9.8*(1000-683)/0.071);
eq2=sin(0.959931-x)*besselk(0,(9.8*(1000-683)/0.071)^ ...
    0.5*0.000115*sin(x))+(9.8*(1000-683)/0.071)^0.5*y* ...
    besselk(1,(9.8*(1000-683)/0.071)^0.5*0.000115*sin(x));
[x,y]=solve(eq1,eq2)
得到的结果为：
x =
7.2408611253328779611417644360562
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-2*2480)/ ...
    (3*(1000-683))-y*(sin(x))^2/0.000115-2*sin(x)* ...
    sin(x-0.959931)/0.000115^2/(9.8*(1000-683)/0.071);
f(2)=sin(0.959931-x)*besselk(0,(9.8*(1000-683)/0.071)^ ...
    0.5*0.000115*sin(x))+(9.8*(1000-683)/0.071)^0.5*y* ...
    besselk(1,(9.8*(1000-683)/0.071)^0.5*0.000115*sin(x));
得到的结果：
x = 0.957676, y = -8.58706e-007
Eq1(x, y) = -5.15143e-014, Eq1(x, y) = -1.38778e-017
我觉得没有太多需要的说明的,所以就不多写了,
 
另外,我注意到,前面用符号数学工具箱求出的解与使用fsolve得到的结果相差2*pi,有兴趣可自行验证.
 
由于百度知道系统抽风,另一个问题的答案“正在提交中”,我把主要内容贴到这里,供参考.
 
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个函数完全一致,但支持微分、积分等符号运算.