求助matlab大神，怎样把一个变量带入3元方程中，做出3维图形~~~~~~~~~~~~~~~~

近段时间不知道怎么原因，校园网不能访问百度空间，所以你的问题我没法看。
后来用手机打开看的，虽然费了一些波折，但总还是看到了。

对这个问题，其实用我之前在下面这个提问中说到的方法解决并不困难：
http://zhidao.baidu.com/question/561545056?&oldq=1
当时我说到“无解时会报错（可通过try...catch结构处理）”，但没进一步做，现在正好通过这个例子给你演示一下。

顺便问一下，现在的f(x,y,z)=0方程与之前在这个提问中的不一样？
http://zhidao.baidu.com/question/561281304?&oldq=1
我开始的时候用上面说的方法去做，求出来的结果画出图来看上去和以前的差别很大，还以为是方法有问题呢，后来仔细看，原来方程变了。

仍然用这个帖子里的方法：
http://zhidao.baidu.com/question/561151606?&oldq=1
取不同的z值画x-y曲线，现在的方程f(x,y,z)=0得到的一组x-y曲线如下：

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

下面是程序和绘图结果（全部代码保存一个文件中运行）。
function zdN= 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;clfmesh(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 zf=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结构避免方程无解时出错，如果无解，则返回NaNtry 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

楼上的程序写的还不错，但是，建议在提交答案之前至少应该做一些基本的验证工作。毕竟，提供一个错误的答案比起没有答案不见得会好到哪里去。

我对这个问题的直观反应是，两个曲面的交线应该会比较平滑，而不应该像楼上给出来的这样杂乱无章。不妨通过一小段绘图的代码验证一下：
N = 25;x=linspace(0.52, 0.53, N); y=linspace(-1e-9, 0, N); z=linspace(1, 20000, N);[x,y,z]=meshgrid(x,y,z);F1=-2.*pi.*0.05415.*0.0000002.*sin(x).*sin((5.*pi./6)+x)-4./3.*pi.* ... 0.0000002.^3.*2000.*z.*9.78+2./3.*pi.*z.*9.78.*0.0000002.^3.* ... (1820-1000).*(cos(x)).^3-2./3.*pi.*z.*9.78.*0.0000002.^3.* ... (1820+1000)-pi.*z.*9.78.*0.0000002.^2.*(y+0.0000002.*cos(x)).* ... (1820-1000).*(sin(x)).^2; F2=-sin((5.*pi./6)+x).*besselk(0,(z.*9.78.*(1820-1000)./0.05415).^ ... 0.5.*0.0000002.*sin(x))+(z.*9.78.*(1820-1000)./0.05415).^0.5.* ... y.*besselk(1,(z.*9.78.*(1820-1000)./0.05415).^0.5.*0.0000002.*sin(x));clf; p = patch(isosurface(x,y,z,F1,0,z));set(p,'facecolor','w','EdgeColor','flat');p = patch(isosurface(x,y,z,F2,0,z));set(p,'facecolor','interp','EdgeColor','interp');view(-10,35);grid on;xlabel('x')ylabel('y')zlabel('z')绘制得到的曲面图形如下。



从图中我们很直观地初步得到几点结论：

（1）方程1的曲面对y几乎不相关。这一点其实比较容易理解，因为方程中y只出现在下面这一项中
(y+0.0000002.*cos(x))我们把这个式子的两项绝对值做一下对比：
max(abs(y(:)./(0.0000002.*cos(x(:)))))可以得到，最大值仅为0.0058，所以y在这个方程中完全可以忽略。

（2）方程1的x和z之间几乎是线性关系。

（3）两个曲面的交线近似是一条很平滑的空间曲线。

接下来，我们尝试画出这条曲线。代码如下：
syms x y zeq1=-2.*pi.*0.05415.*0.0000002.*sin(x).*sin((5.*pi./6)+x)-4./3.*pi.* ... 0.0000002.^3.*2000.*z.*9.78+2./3.*pi.*z.*9.78.*0.0000002.^3.* ... (1820-1000).*(cos(x)).^3-2./3.*pi.*z.*9.78.*0.0000002.^3.* ... (1820+1000)-pi.*z.*9.78.*0.0000002.^2.*(y+0.0000002.*cos(x)).* ... (1820-1000).*(sin(x)).^2; eq2=-sin((5.*pi./6)+x).*besselk(0,(z.*9.78.*(1820-1000)./0.05415).^ ... 0.5.*0.0000002.*sin(x))+(z.*9.78.*(1820-1000)./0.05415).^0.5.* ... y.*besselk(1,(z.*9.78.*(1820-1000)./0.05415).^0.5.*0.0000002.*sin(x));Z=solve(eq1,z);eq = subs(eq2, z, Z);hold onh = ezplot(eq, [0.52 0.53], [-1e-9 0]);X = get(h,'XData');Y = get(h,'YData');view(-10,35)grid onz = subs(Z, {x y}, {X Y});x(z>20000)=NaN;y(z>20000)=NaN;z(z>20000)=NaN;set(h, 'XData', X, 'YData', Y, 'ZData', z, 'linewidth', 2);title('')p=findall(gca,'type','patch');a=0.2;set(p,'Facealpha',a,'EdgeAlpha',a);得到的绘图结果如下：

对上面的代码，简单说明几点：
（1）可以单独运行，也可以接在前面的代码后面运行（单独运行时不画曲面）；
（2）
（3）为了突出显示交线，把两个曲面设为半透明，如果不想要这两个曲面，可以直接在把代码的最后两行换成：
delete(p)axis tight即可删除曲面而只保留曲线。

以上程序在2007b上测试通过，其它环境暂未测试。
全部代码已作为附件上传，希望对楼主有帮助，如有问题可追问。

其实这个问题比较简单：在给定区间范围内绘制v=f(z,n)的曲面，公式是现成的，用meshgrid生成区域的座标点，然后代入公式计算，最后用mesh或surf之类的函数绘图即可。

唯一的问题是，v不仅是z和n的函数，还有一个自变量x。但x并非真正的独立变量，而是受另外两个三元方程约束的，对于给定的z，可以用question/577643560.html中的方法求出对应的x来。

你在空间贴出的代码是上述问题中的第一种方法。如果是用Maple内核的MATLAB（例如6.5或2007b，这两个版本上我做了测试），应考虑使用第二种方法，精度更高。

以下是代码，供参考（已作为附件上传）：

% 绘图座标范围
ni = linspace(2, 100, 20);
zi = linspace(1, 20000, 30);

% 对于给定的z，求解对应的x
syms x y z 
eq1=-2.*pi.*0.05415.*0.0000002.*sin(x).*sin((5.*pi./6)+x)-4./3.*pi.* ... 
    0.0000002.^3.*2000.*z.*9.8+2./3.*pi.*z.*9.8.*0.0000002.^3.* ... 
    (1820-1000).*(cos(x)).^3-2./3.*pi.*z.*9.8.*0.0000002.^3.* ... 
    (1820+1000)-pi.*z.*9.8.*0.0000002.^2.*(y+0.0000002.*cos(x)).* ... 
    (1820-1000).*(sin(x)).^2;  
eq2=-sin((5.*pi./6)+x).*besselk(0,(z.*9.8.*(1820-1000)./0.05415).^ ... 
    0.5.*0.0000002.*sin(x))+(z.*9.8.*(1820-1000)./0.05415).^0.5.* ... 
    y.*besselk(1,(z.*9.8.*(1820-1000)./0.05415).^0.5.*0.0000002.*sin(x)); 
xi = sym( zi*0 ); 
yi = xi; 
for i = 1 : length(zi) 
    z0 = zi(i); 
    [xi(i), yi(i)] = solve(subs(eq1,z,z0), subs(eq2,z,z0)); 
end

% SOLVE的求解结果x位于要求的坐标范围之外，将其减小2*pi并转为数值类型
xi = double(xi-2*pi);

% 计算函数值
[n, z] = meshgrid(ni, zi);
[n, x] = meshgrid(ni, xi);
v = 2*0.05415*(-tan(5*pi/6+x-pi)./(((z*9.8*(1820-1000)/0.05415) ...
    .^0.5).*besselk(1,((z*9.8*(1820-1000)/0.05415).^0.5)*0.0000002 ...
    .*sin(x)))).^2.*((z*9.8*(1820-1000)/0.05415).^0.5).*besselk(1, ...
    ((z*9.8*(1820-1000)/0.05415).^0.5)*0.0000002.*n).*(1-1/3.* ...
    n.^-1+n.^-3-15/4.*n.^-4-4.46/1000.*(n-1.7).^-2.867)./(3*0.8937 ...
    /1000*0.0000002*1.0);


% 绘制曲面及标注图形
mesh(n, z, v)
view(30,15)
xlabel('n')
ylabel('z')
zlabel('v')

曲面效果如下图所示。