matlab中怎样从曲线中获得精确的坐标值~~~~~~~~~~~~~~
提供几种不同的做法，供参考。
方法1：
=====
直接从绘图数据插值（经检验z数据是单调增加的），代码如下：
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.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));
Z=solve(eq1,z);
eq=subs(eq2,z,Z);
h=ezplot(eq,[0.52?0.5245],[-1e-9?0]);
X=get(h,'XData');
Y=get(h,'YData');
view(-10,35)
grid?on
zz=subs(Z,{x?y},{X?Y});
X(zz>20000)=[];
Y(zz>20000)=[];
zz(zz>20000)=[];
set(h,?'XData',?X,?'YData',?Y,?'ZData',?zz,?'linewidth',?2);
title('')
zlabel('z')
axis?tight
box?off
zi?=?[1?10?100?1000?5000?10000?19000?20000];
xi?=?interp1(zz,?X,?zi,?'s');
yi?=?interp1(zz,?Y,?zi,?'s');
vpa([xi;?yi].',10)
举例：
------
上面的代码中z分别取[1 10 100 1000 5000 10000 19000 20000]，得到的结果如下（每行代表一组对应的x、y）：
[?.5235988072,?.1770687259e-12][?.5235990919,?-.3682753278e-13]
[?.5236019385,?-.2134703182e-11]
[?.5236304028,?-.2157265249e-10]
[?.5237568792,?-.9522223488e-10]
[?.5239149030,?-.1794854554e-9]
[?.5241991463,?-.3217413726e-9]
[?.5242307131,?-.3370566367e-9]
把数据代回方程，看一下精度：
>>?eq?=?subs([eq1;?eq2],?{x?y?z},?{xi?yi?zi})eq?=
1.0e-005?*0.0000-0.0000-0.0000-0.00000.0000-0.0000-0.0000-0.0000
0.20960.25190.38850.01820.00550.00720.0034-0.0010
>>?norm(eq)
ans?=
5.0868e-006说明：
------
（1）你所贴代码的前半部分本来用于说明两个曲面相交的，如果只画线，该部分属于多余；
（2）顺便修正一点小错误——原代码的这两句有点小问题：
x(z>20000)=NaN;y(z>20000)=NaN;
我上次还有点奇怪，为什么坐标范围看上去不太对劲，但没有深究。刚才发现是这两句写错了，x、y应该是大写的。另外，现在为符合插值的需要应把右侧NaN改为空矩阵（[]）。
（3）有多种插值方法可用，我使用了最有利于平滑曲线的样条插值，楼主也可以试试其它做法。
（4）注意z0取值的范围最好不要超过[min(z)?max(z)]=[8.45?19943]，如果超出该范围，部分方法需要通过参数指定允许interp1外插（默认情况下，样条插值允许外插，但线性插值不允许）。
（5）由于方法本身是基于部分离散点信息得到未知点的，所以很难说是否能保证精度。
（6）兼容性：这里的绘图代码在2007b上测试没问题，但在6.5上不行，更高版本没试，也可能会有问题。
方法2：
=====
从原始问题出发，直接解方程，代码如下：
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.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));
zi?=?[1?10?100?1000?5000?10000?19000?20000];
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
vpa([xi-2*pi;?yi].',10)
举例：
------
与方法1 的测试用例相同，得到的结果如下：
[.523598806,?-.3251499235e-13][.523599091,?-.2887359289e-12]
[.523601938,?-.2523219710e-11]
[.523630402,?-.2159081181e-10]
[.523756878,?-.9522774887e-10]
[.523914902,?-.1794926954e-9]
[.524199145,?-.3217447875e-9]
[.524230712,?-.3370556700e-9]
同样，检验一下精度：
>>?eq?=?subs([eq1;?eq2],?{x?y?z},?{xi-2*pi?yi?zi});>>?eq=double(eq)
eq?=
1.0e-018?*0.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00020.00210.01820.08120.15390.27740.2907
>>?norm(eq)
ans?=
4.3823e-019说明：
------
（1）由于求解方程组得到x总是位于要求的坐标范围之外，所以把它减小2*pi。
（2）这种方法需要符号数学工具箱的支持，精度远高于前一种方法。
（3）兼容性：用solve解方程的结果可能和符号数学工具箱版本相关，我用Maple内核的6.5和2007b做了测试，可以求出想要的结果，但在MuPad内核的版本上可能会有问题。
方法3：
=====
使用数值方法解方程。本来以为第2种方法有些条件下会失效，所以考虑这种做法，但从实际情况看，方法2的效果很好，所以，这个就不做了。
====================================
最后，八卦几句：
已经多次回答楼主的问题了，有点好奇您是从哪里找到这么多复杂的表达式要求解的，和您的工作有关吗？而且难得的是，问题看上去虽然比较复杂，但根据楼主提供的信息，刚好都是可以求解的，简直就像是考试的试题一般。
顺便说一下，我注意到，这次的方程与上次我回答的问题相比，有7处由9.78换成了9.8（想来应该是重力加速度），应该是有意为之的吧？