近段时间不知道怎么原因,校园网不能访问百度空间,所以你的问题我没法看.

后来用手机打开看的,虽然费了一些波折,但总还是看到了.

对这个问题,其实用我之前在下面这个提问中说到的方法解决并不困难：

当时我说到“无解时会报错(可通过try...catch结构处理)”,但没进一步做,现在正好通过这个例子给你演示一下.

顺便问一下,现在的f(x,y,z)=0方程与之前在这个提问中的不一样?

我开始的时候用上面说的方法去做,求出来的结果画出图来看上去和以前的差别很大,还以为是方法有问题呢,后来仔细看,原来方程变了.

仍然用这个帖子里的方法：

取不同的z值画x-y曲线,现在的方程f(x,y,z)=0得到的一组x-y曲线如下：

从图中可见,x和y之间有很好的线性程度(甚至斜率都差别不大).这个特点也可以作为你进一步研究的参考.

下面是程序和绘图结果(全部代码保存一个文件中运行).

function zd

N= 10;

x = linspace(eps, 0.000001-eps, N);

z = linspace(1,50000, N);

[x, z] = meshgrid(x, z);

y = arrayfun(@f_xz, x, z);

% 方程1：R=f(x,y,z)

R=(-2.*pi.*0.05415./(2.*x)).*(1./2.*(((1820-1000).*z.*9.8./...

0.05415)).*(y.^2.*(x-y./3)-4./3.*((2000-1000)./(1820-...

1000)).*x.^3-(((2.*x-y).*y).^0.5).^2.*((-2.* ((((2.*x-...

y).*y).^0.5)./x).*(((((2.*x-y).*y).^0.5)./x.*cos(5.*pi...

/6)+((x.^2-(((2.*x-y).*y).^0.5).^2).^0.5./x).*sin(5.*pi...

/6))).*0.05415./(1820.*x.^2.*z.*9.8)-(4./3).*2000./1820+...

(2./3-((x.^2-(((2.*x-y).*y).^0.5).^2).^0.5./x)+(1./3).*...

(((x.^2-(((2.*x-y).*y).^0.5).^2).^0.5./x)).^3)+1000./...

1820.*(2./3+((x.^2-(((2.*x-y).*y).^0.5).^2).^0.5./x)-...

(1./3).*(((x.^2-(((2.*x-y).*y).^0.5).^2).^0.5./x)).^3))...

.*x./((1000/1820-1).* ((((2*x-y).*y).^0.5)./x).^2)))).^2;

clf

mesh(x, z, R);

view(-150,35)

xlabel('x')

ylabel('z')

zlabel('R')

% 由方程2：f(x,y,z)=0,计算y=g(x,z)

function y0 = f_xz(x0, z0)

syms x y z

f=1/2*((1820-1000)*z*9.8/0.05415)*(y^2*(x-y/3)-4/3*((2000-...

1000)/(1820-1000))*x^3-(((2*x-y)*y)^0.5)^2*((-2*((((2*...

x-y)*y)^0.5)/x)*((((2*x-y)*y)^0.5)/x*cos(5*pi/6)+((x^2-...

(((2*x-y)*y)^0.5)^2)^0.5/x)*sin(5*pi/6))*0.05415/(1820*...

x^2*z*9.8)-(4/3)*2000/1820+(2/3-((x^2-(((2*x-y)*y)^0.5)...

^2)^0.5/x)+(1/3)*(((x^2-(((2*x-y)*y)^0.5)^2)^0.5/x))^3)...

+1000/1820*(2/3+((x^2-(((2*x-y)*y)^0.5)^2)^0.5/x)-(1/3)...

*(((x^2-(((2*x-y)*y)^0.5)^2)^0.5/x))^3))*x/((1000/1820-...

1)* (((((2*x-y)*y)^0.5)/x))^2)))-(((2*x-y)*y)^0.5*...

sin(acos(y/x-1)-(5*pi/6)));

% 对给定x、z,求解y

% 使用try...catch结构避免方程无解时出错,如果无解,则返回NaN

try

y_xz = @(x0,z0) fzero( @(y0) double(subs(f,{x,y,z}, ...

{x0,y0,z0})), [0 x0], optimset('TolFun', eps));

y0 = y_xz(x0, z0);

catch

disp(lasterr);

y0 = NaN;

end