管道界面面积分割问题 （在线等）

gy8409T

一个通风管道，界面是个半径为1的圆，用隔板分成上下两部分，上部分高为h，给定误差0.001米，求h的数值解，是上下两部分截面面积比为1：2（写出算法的程序和框图）

明心

<>运算和函数符号采用matlab 标准</P><>利用几何或定积分，得到上部分截面的面积s=(acos(1-h)-(1-h).*(1-(1-h).^2).^0.5)</P><>圆面积=pi</P><P>1/3*圆面积＝上部分截面的面积</P><P>求h高度转化为求方程y=1/3*pi-(acos(1-h)-(1-h).*(1-(1-h).^2).^0.5)的零点</P><P>1：构造函数m文件f1_.m</P><P>**************************************************************************</P><P>function y=f1_(h)
y=1/3*pi-(acos(1-h)-(1-h).*(1-(1-h).^2).^0.5);</P><P>***************************************************************************</P><P>2：利用作图，找到方程的大致零点，matlab程序是：</P><P>******************************************************************************************</P><P>h=0:0.01:1;
&gt;&gt; plot(h,f1_(h))
&gt;&gt; hold on
&gt;&gt; plot(h,zeros(length(h),1))</P><P>*******************************************************************************************</P><P>发现零点h在0.7~0.8之间</P><P>3：利用二分法求零点</P><P>相应的m文件是</P><P>********************************************************************************************</P><P>% bisec_n(f_name, a,c): bisection method without graphics
%      f_name: definition of the equation to solve
%       a, c : end points of initial interval
% Copyright S. Nakamura, 1995
function bisec_n(f_name, a,c)
f_name
%       tolerance : tolerance
%       it_limit : limit of iteration number
%       Y_a, Y_c : y values of the current end points
fprintf( 'Bisection Scheme\n\n' );
tolerance = 0.000001;  it_limit = 30;
fprintf( ' It.  a           b           c            f(a)   ');
fprintf( '     f(b)        f(c)\n' );
it = 0;
Y_a = feval(f_name, a ); Y_c = feval(f_name, c );
if ( Y_a*Y_c &gt; 0 )  
     fprintf( '\n\n Stopped because   f(a)f(c) &gt; 0 \n' );
     fprintf( '\n Change a or b and run again.\n' );
else   
   while 1
      it = it + 1;
      b = (a + c)/2;  Y_b = feval(f_name, b );
      fprintf('%3.0f %10.6f, %10.6f', it, a, b );
      fprintf('%10.6f,  %10.6f, %10.6f, %10.6f\n', c, Y_a, Y_b, Y_c );
      if ( abs(c-a)&lt;=tolerance )
         fprintf( '  Tolerance is satisfied. \n' );break
      end
      if ( it&gt;it_limit )
        fprintf( 'Iteration limit exceeded.\n' ); break
      end  
      if( Y_a*Y_b &lt;= 0 )     c = b;    Y_c = Y_b;
      else                   a = b;    Y_a = Y_b;
      end
   end
   fprintf('Final result: Root = %12.6f \n', b );
end</P><P>**************************************************************************</P><P>调用格式</P><P>*************************************************************************</P><P>bisec_n('f1_',0.7,0.8)</P><P>**************************************************************************</P><P>4:结果</P><P>h=0.735068
</P>

gy8409T

[em01]谢谢楼上