问题我已经贴到百度空间了,请大神帮我看看,来日一定报答

问题我已经贴到百度空间了,请大神帮我看看,来日一定报答
问题我已经贴到百度空间了,请大神帮我看看,来日一定报答

问题我已经贴到百度空间了,请大神帮我看看,来日一定报答

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

 

我对这个问题的直观反应是,两个曲面的交线应该会比较平滑,而不应该像楼上给出来的这样杂乱无章.不妨通过一小段绘图的代码验证一下：

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 z
eq1=-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 on
h = ezplot(eq, [0.52 0.53], [-1e-9 0]);
X = get(h,'XData');
Y = get(h,'YData');
view(-10,35)
grid on
z = 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上测试通过,其它环境暂未测试.

全部代码已作为附件上传,希望对楼主有帮助,如有问题可追问.