matlab实现曲面积分（surf_integral函数）

深碍√TFBOYSˉ_ 2022-11-25 00:40 199阅读 0赞

MATLAB语言并未直接提供曲面积分的现成函数，因此，此处给出计算曲面积分的函数。

### 目录 ###

*   *   *  函数说明
         *  应用举例
         *   *  第一类曲面积分
             *  第二类曲面积分
         *  函数实现

### 函数说明 ###

function I = surf_integral(f,vars,t,a,b)
    %surf_integral
    %第一类曲面积分
    %   I = surf_integral(f, z, [x,y], [x_m,x_M], [y_m,y_M])
    %   I = surf_integral(f, [x,y,z], [u,v], [u_m,u_M], [v_m,v_M])
    %  Examples:
    %  计算int_int(x^2*y+z*y^2)dS, 积分曲面如下：
    %  x=ucosv, y=usinv, z=v, 0<=u<=a, 0<=v<=2*pi
    %  MATLAB求解语句
    %  syms u v; syms a positive;
    %  x=u*cos(v); y=u*sin(v); z=v; f=x^2*y+z*y^2;
    %  I = surf_integral(f,[x,y,z],[u,v],[0,a],[0,2*pi])
    %
    %第二类曲面积分
    %   I = surf_integral([P,Q,R], z, [x,y], [x_m,x_M], [y_m,y_M])
    %   I = surf_integral([P,Q,R], [x,y,z], [u,v], [u_m,u_M], [v_m,v_M])
    %  注意：I = int_int_S(P*dydz+Q*dxdz+R*dxdy)
    %  Examples:
    %  计算曲面积分int_int(x*y+z)dxdy, 积分曲面如下：
    %  (x/a)^2+(y/b)^2+(z/c)^2=1的上半部，且积分沿椭球面的上面。
    %  引入参数方程：x=a*sin(u)*cos(v)，y=b*sin(u)*sin(v)，z=c*cos(u);
    %  且0<=u<=pi/2, 0<=v<=2*pi 。
    %  MATLAB求解语句
    %  syms u v; syms a b c positive;
    %  x=a*sin(u)*cos(v); y=b*sin(u)*sin(v); z=c*cos(u);
    %  I = surf_integral([0, 0, x*y+z],[x,y,z],[u,v],[0,pi/2],[0,2*pi])

### 应用举例 ###

#### 第一类曲面积分 ####

1.  计算  ∬ ( x 2 y + z y 2 ) d S \\iint(x^2y+zy^2)dS ∬(x2y\+zy2)dS, 积分曲面如下：  
     x = u c o s v ,    y = u s i n v ,    z = v ,    0 ≤ u ≤ a ,    0 ≤ v ≤ 2 π x=ucosv,~~ y=usinv, ~~z=v, ~~0\\leq u\\leq a,~~ 0\\leq v\\leq 2\\pi x=ucosv,  y=usinv,  z=v,  0≤u≤a,  0≤v≤2π

**求解方法**

syms u v; syms a positive;
    x=u*cos(v); y=u*sin(v); z=v; f=x^2*y+z*y^2;
    I = surf_integral(f,[x,y,z],[u,v],[0,a],[0,2*pi])

#### 第二类曲面积分 ####

1.  计算曲面积分  ∬ ( x y + z ) d x d y \\iint(xy+z)dxdy ∬(xy\+z)dxdy, 积分曲面如下：  
     ( x a ) 2 + ( y b ) 2 + ( z c ) 2 = 1 (\\frac\{x\}\{a\})^2+(\\frac\{y\}\{b\})^2+(\\frac\{z\}\{c\})^2=1 (ax)2\+(by)2\+(cz)2=1 的上半部，且积分沿椭球面的上面。  
    注：引入参数方程： x = a s i n ( u ) c o s ( v ) ， y = b s i n ( u ) s i n ( v ) ， z = c c o s ( u ) ; x=a sin(u) cos(v)，y=b sin(u) sin(v)，z=ccos(u); x=asin(u)cos(v)，y=bsin(u)sin(v)，z=ccos(u); 且  0 ≤ u ≤ π 2 , 0 ≤ v ≤ 2 π 0\\leq u \\leq \\frac\{\\pi\}\{2\}, 0\\leq v\\leq 2\\pi 0≤u≤2π,0≤v≤2π

**求解方法**

syms u v; syms a b c positive;
    x=a*sin(u)*cos(v); y=b*sin(u)*sin(v); z=c*cos(u);
    I = surf_integral([0, 0, x*y+z],[x,y,z],[u,v],[0,pi/2],[0,2*pi])

### 函数实现 ###

function I = surf_integral(f,vars,t,a,b)
    if length(f)==1
        if length(vars)~=1
            E = simplify(sum(diff(vars,t(1)).^2));
            F = sum(diff(vars,t(1)).*diff(vars,t(2)));
            G = simplify(sum(diff(vars,t(2)).^2));
        else
            E = simplify(1+diff(vars,t(1))^2);
            F = diff(vars,t(1))*diff(vars,t(2));
            G = simplify(1+diff(vars,t(2))^2);
        end
        I = int(int(simplify(f*sqrt(E*G-F^2)),t(1),a(1),a(2)),t(2),b(1),b(2));
    else
        if length(vars)~=1
            A = diff(vars(2),t(1))*diff(vars(3),t(2)) - diff(vars(3),t(1))*diff(vars(2),t(2));
            B = diff(vars(3),t(1))*diff(vars(1),t(2)) - diff(vars(1),t(1))*diff(vars(3),t(2));
            C = diff(vars(1),t(1))*diff(vars(2),t(2)) - diff(vars(2),t(1))*diff(vars(1),t(2));
        else
            A = - diff(vars,t(1));
            B = - diff(vars,t(2));
            C = 1;
        end
        f = f(:); abc = [A, B, C];
        I = int(int(simplify(abc*f),t(1),a(1),a(2)),t(2),b(1),b(2));
    end

--------------------

此函数源文件可前往下面网址下载：

[surf\_integral.m下载通道][surf_integral.m]

[surf_integral.m]: https://download.csdn.net/download/weixin_43964993/12692475