int .

“msharpmath”, The Simple is the Best

[100] 003 Tutorial int                     integration

Click the above button for full-text tutorial of   ‘int’ .

//==================================================================================
//   Umbrella ‘int’   for integration
//==================================================================================

//==================================================================================
//   Single Integration
//==================================================================================
//    int .x[ n = 26, g = 1 ] (a,b) ( <<opt>>, f(x), g(x), h(x), … )
//    x = (a,b).span(n,g = 1);  int [x] ( <<opt>>, f(x), g(x), h(x), … )

n = number of points used in integration
g = ratio between two succesive grid size, dx_i = g dx_(i-1), normally g = 1 by default
a = start point of integration
b = end point of integration, a != b but it is not necessary to be a < b.
f(x), g(x), h(x) = integrand for integration

//———————————————————————————-
//  .poly(nth)        Newton-Cotes integration with polynomial of degree nth <= 8
//  .quad(nth)       Gauss-Legendre quadrature with polynomial of degree nth <= 5
//  .togo(F)             takeout an integrated function as a matrix F
//  .return               return an integrated function as a matrix
//  .plot                   plot integrated function for one-variable case
//———————————————————————————-

#> int .x[5](0,1) ( 1/(1+x*x) ); pi/4;       // a default method of integration
ans =      0.78539816
ans =      0.78539816

//————————————
// Newton-Cotes method,                .poly(nth)  nth is the degree of polynomial
//————————————
// n = 5 points at x = 0, 0.25, 0.5, 0.75, 1 (a = 0, b = 1)
// note that 4 steps in 5 points

#> int .x[5](0,1) ( exp(2*x) ).poly(0);
#> 0.25*( exp(2*0) + exp(2*0.25) + exp(2*0.5) + exp(2*0.75) ) ;
ans =        2.462173
ans =        2.462173

#> int .x[5](1,0) ( exp(2*x) ).poly(0);        // reverse direction from 1 to 0
#> -0.25*( exp(2*1) + exp(2*0.75) + exp(2*0.5) + exp(2*0.25) );
ans =      -4.0594371
ans =      -4.0594371

#> int .x[5](0,1) ( exp(2*x) ).poly(1);        // trepezoidal rule with a line
#> 0.25/2*( exp(2*0) + 2*exp(2*0.25) + 2*exp(2*0.5) + 2*exp(2*0.75) + exp(2*1) ) ;
ans =       3.2608051
ans =       3.2608051

#> int .x[3](0,1) ( exp(2*x) ).poly(2);        // Simpson’s 1/3 rule
#> 0.25/3*( (exp(2*0)+4*exp(2*0.25)+exp(2*0.5)) + (exp(2*0.5)+4*exp(2*0.75)+exp(2*1)) ) ;
ans =       3.1956051
ans =       3.1956051

//————————————
// Gauss quadrature
//————————————
#> int .x[2](-1,1) ( exp(x) ).quad(1);     2*exp(0);
ans =               2
ans =               2

#> int .x[2](-1,1) ( exp(x) ).quad(2);      exp(-1/sqrt(3)) + exp(1/sqrt(3));
ans =       2.3426961
ans =       2.3426961

#> int .x[2](-1,1) ( exp(x) ).quad(3);     5/9*( exp(-sqrt(0.6))+exp(sqrt(0.6)) )+8/9*exp(0);
ans =       2.3503369
ans =       2.3503369

//==================================================================================
//  Integrated Function and Plot
//==================================================================================
//  int .x[n=26,g=1](a,b) ( <<opt>>, f(x), g(x), … ) .togo(F,G, …);
//  int .x[n=26,g=1](a,b) ( <<opt>>, f(x), g(x), … ) .return;

#> int .x[5](0,1) ( x, x*x ) .togo(P,Q);  P;  Q;
ans =  [           0.5      0.333333 ]
P =
[             0 ]
[       0.03125 ]
[         0.125 ]
[       0.28125 ]
[           0.5 ]
Q =
[             0 ]
[    0.00520833 ]
[     0.0416667 ]
[      0.140625 ]
[      0.333333 ]

#> int .x[5](0,1) ( x, x*x ) .return;
ans =
[             0             0 ]
[       0.03125    0.00520833 ]
[         0.125     0.0416667 ]
[       0.28125      0.140625 ]
[           0.5      0.333333 ]

//———————————–
// symmetric plot
//———————————–
#> x = (0,8*pi) .span(1001);;
#> int [x] ( sin(x^1.5), cos(x^1.5) ) .togo( S, C );;
#> plot( [ -S.rev, S ], [ -C.rev, C ] );

//==================================================================================
//  Multiple Integration
//==================================================================================
//    int .x[n=11,g=1](a,b) .y[n=11,g=1](c(x),d(x)) ( <<opt>>, f(x,y) )
//    int .x[n= 6,g=1](a,b) .y[n= 6,g=1](c(x),d(x)) .z[n= 6,g=1](u(x,y),v(x,y)) ( <<opt>>, f(x,y,z) )


//———————————————————————————-

#> int .x[11](0,1) .y[11]( -sqrt(x-x*x),sqrt(x-x*x) ) ( x-x*x-y*y );   pi/32;
ans =     0.098174621
ans =      0.09817477

#> int .x[51](0,1) .y[51]( -sqrt(x-x*x),sqrt(x-x*x) ) ( x-x*x-y*y ).poly(1);
ans =     0.098131619

#>   int .x(0,3).z(1,x).y(z-x,z+x) ( exp(2*x)*(2*y-z) );
#>   1/8 + 59*exp(6)/8;
ans =       2975.4124
ans =       2975.4124

#>   int .x[26](0,3) .z[26](1,x) .y[26](z-x,z+x) ( exp(2*x)*(2*y-z) ).poly(1);  // trepezoidal rule
#>   1/8 + 59*exp(6)/8;
ans =       3011.1664
ans =       2975.4124

As expected, the Gauss quadrature is more accurate than the trepezoidal rule, when the integrand function is given analytically.

//==================================================================================
//  Area Properties for Closed Curves in 2D
//==================================================================================
//   .int1da(a)     area enclosed by a closed curve
//   .intxda(ax)    area properties
//   .intyda(ay)
//   .intxxda(axx)
//   .intxyda(axy)
//   .intyyda(ayy)


//———————————————————————————-

%> example 3-4-1  area and moment of a closed curve ========================
#> .hold;
#>   plot .line(-2,0, 2,0);
#>   plot .x(2,-2) ( 4-x*x ) .link .int1da(area) .intyyda(ayy);
#>   plot .line(0,-1, 0,4);
#> plot;
#> area / (32/3);      ayy / ( 4096/105 );    // ratio with an exact value
ans =      0.99960011
ans =      0.99906713


#> example 3-4-2  area and moment of a closed curve =========================
#> .hold;
#>   plot .line(-2,0, 2,0);
#>   plot .x(2,-2) ( 4-x*x ) .link .zrot(pi/6) .int1da(area) .intyyda(ayy);
#>   plot .line(0,-1, 0,4);
#>   plot .line(0,0, 0,4) .zrot(pi/6);
#> plot;
#> area / (32/3); ayy / ( 4096/105 );     // after rotation, a_yy has changed
ans =       0.9996002
ans =       0.8039512

%> first and second moment ========================================
plot .ellipse[101](3,2)   // 101 points are used (100 equal spacing)
.int1da(a).intxda(ax).intyda(ay)
.intxxda(axx).intxyda(axy).intyyda(ayy);

#> a/(6*pi); ax; ay; axx; axy; ayy;  (axx+ayy) / (39*pi/2);  // 0.13% error
ans =      0.99934218
ax =  6.1552078e-007
ay = -8.7607567e-008
axx =       42.355721
axy = -9.4752596e-007
ayy =        18.82476
ans =      0.99868471

//==================================================================================
//  Surface Integral
//==================================================================================
// .area(a)
// .intda (f(x,y,z), …) (fa, …)


//———————————————————————————-
// .intdA * (Fx(x,y,z),Fy(x,y,z),Fz(x,y,z)) (fscalar)
// .intdA ^ (Fx(x,y,z),Fy(x,y,z),Fz(x,y,z)) (Fmatrix)
// .intdA   (f(x,y,z), …)                 (Fmatrix, …)


//———————————————————————————-

%> example 3-5-1 =================================================

#>  plot .x[51](-2,2) ( 0, 4-x*x ) .fill[51]
.area(a)
.intda ( x, y, x*x, x*y, y*y, x*x*(1+0.1*(x*y)) ) ( ax,ay,axx,axy,ayy, rhoaxx );

#>  a; ax; ay; axx; axy; ayy; rhoaxx;
a =         10.6624
ax =  9.0205621e-017
ay =       17.049604
axx =        8.533332
axy =  7.0513799e-017
ayy =       38.951029
rhoaxx =        8.533332

%> example 3-5-2 =================================================
#> plot .x[31](1,sqrt(5)) (0) .zrev[361](0,2*pi)
.area(a)
.intda( x*x+y )(b);

#> a/(4*pi);  b/(6*pi);
ans =       0.9999493
ans =      0.99973176

%> example 3-5-3  intdA * ( Fx,Fy,Fz )(f)  =================================

#> plot .sphere[101,101](1) .intdA*(x,2*y,3*z)(fs);

#> fs/(8*pi);
ans =      0.99911615

 

%> example 3-5-4  intdA * ( Fx,Fy,Fz )(f)  =================================

#> plot .sphere[101,101](1) .intdA*(x,y,z)(fs);

#>  fs/(4*pi);
ans =       0.9990956

 

%> example 3-5-5  intdA ^ ( Fx,Fy,Fz )(Fmatrix)  ============================

#> plot .disk[101,101](1) .intdA^(sqrt(x*x+y*y),x*x+y*y,0)(F);

#> F;  [-pi/2, 2*pi/3, 0];

F =  [      -1.56814       2.09193             0 ]
ans =  [       -1.5708        2.0944             0 ]

 

%> example 3-5-6  intdA  ( f )(Fmatrix)   ==================================

plot .sphere[101,101](1) .intdA (z)(F);

F;  [0,0, 4*pi/3];
F =  [  2.61564e-018  3.32503e-018       4.18535 ]
ans =  [             0             0       4.18879 ]

 

%> example 3-5-7  area and volume of closed surface =========================

#> plot .torus[51,51](0.5,0.5) .area(a) .volume(v);

#> a/pi/pi;  v/(0.25*pi*pi);
ans =      0.99769965
ans =      0.99474528

%> example 3-5-8  area and volume of closed surface =========================

#> plot .frustum[51,5,51](3,2,2) .area(a) .volume(v);

#> a/56; v/(76/3);
ans =               1
ans =       1.0000105

 

%> example 3-5-9  area and volume of closed surface =========================

#> plot .frustum[21,31,21](3,0,4) .area(a) .volume(v);

#> a/(24*pi); v/(12*pi);  // 0.9995395 1.0005861
ans =      0.99489459
ans =      0.99394582

//==================================================================================
//  Line Integral
//==================================================================================
%> intdL*(Fx,Fy,Fz)
%> intdL^(Fx,Fy,Fz), intdL (f)
%> intdl(f,g,…)

==============================================
#>  plot .@t[101](0,2*pi) ( 3*cos(t),2*sin(t),0 )
.intdL*(-y/(x*x+y*y), x/(x*x+y*y) )(fval);

#> fval/(2*pi);
ans =       1.0003291

 

%>  intdL*(Fx,Fy,Fz) ===============================================

#> plot .@t(pi,3*pi) ( cos(t),sin(t),.5*t ) .intdL*(y,x,z) (fval);

fval/(pi*pi);
ans =      0.99999978

 

%> intdL^(Fx,Fy,Fz), intdL (f)  ========================================

#> plot .@t(pi,3*pi) ( cos(t),sin(t),0.5*t ) .intdL^(y,x,z) (Fcrs);

#> plot .@t(pi,3*pi) ( cos(t),sin(t),0.5*t ) .intdL (1) (Fvec);

#> Fcrs; Fvec;
Fcrs =  [  7.72116e-008       3.14159  2.08167e-017 ]
Fvec =  [             0    2.498e-016       3.14159 ]

 

%> intdl(f,g,…)  ==================================================

#> plot .@t[101](pi,3*pi) ( cos(t),sin(t),0.5*t ) .intdl(1, L) (arclen, f);;

#> arclen/(sqrt(5)*pi); f/(2.5*pi*pi);
ans =      0.99986842
ans =       1.0097342

 

Comments are closed.