matlab开发-椭圆积分和函数

matlab开发-椭圆积分和函数。使用AGM算法进行椭圆函数评估。

function [abcn] = agm（a0b0c0tol）
% AGM calculates the Artihmetic Geometric Mean of A and B （see [1]）. 
% The function is used in routines ELLIPJ and ELLIPTIC12.
%
%   [ABCN] = AGM（A0B0C0TOL） carry out the process of the arithmetic geometric 
%   mean starting with a given positive numbers triple （A0B0C0） and returns in 
%   （ABC） the generated sequence. N is a number of steps （returns in the type uint32）.
%
%   The general scheme of the process:
%       A（i） = 1/2*（ A（i-1）+B（i-1） ）;     A（0） = A0;
%       B（i） = sqrt（A（i-1）*B（i-1））;       B（0） = B0;
%       C（i） = 1/2*（ A（i-1）+B（i-1） ）;     C（0） = C0;
%   Stop at the N-th step when A（N） = B（N） i.e. when C（N） = 0.
%
%   See also ELLIPKE ELLIPJ ELLIPTIC12 ELLIPTIC3 THETA.
%
%   References:
%   [1] M. Abramowitz and I.A. Stegun “Handbook of Mathematical
%       Functions“ Dover Publications“ 1965 Ch. 17.1 - 17.6.


if nargin<4 tol = eps; end
if nargin<3 error（‘Not enough input arguments.‘）; end

% pre-allocate space and augment if needed
chunk = 8; mmax = prod（size（a0））;
a = zeros（chunkmmax）;
c = a;
b = a;
a（1:） = a0;
b（1:） = b0;
c（1:） = c0;
n = uint32（ zeros（1mmax） ）;
i = 1;
while any（abs（c（i:）） > tol）
    i = i + 1;
    if i > size（a1）
      a = [a; zeros（chunkmmax）];
      b = [b; zeros（chunkmmax）];
      c = [c; zeros（chunkmmax）];
    end
    a（i:） = 0.5 * （a（i-1:） + b（i-1:））;
    b（i:） = sqrt（a（i-1:） .* b（i-1:））;
    c（i:） = 0.5 * （a（i-1:） - b（i-1:））;
    in = uint32（ find（（abs（c（i:）） <= tol） & （abs（c（i-1:）） > tol）） ）;
    if ~isempty（in）
      [mini] = size（in）;
      n（in） = ones（mini）*（i-1）;
    end
end

% END FUNCTION AGM（）

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        1657  2007-10-23 18:04  agm.m
     文件        3415  2007-10-25 02:17  ellipj.m
     文件        2439  2009-06-20 17:33  ellipji.m
     文件        4888  2009-06-20 17:34  elliptic12.m
     文件        3727  2009-06-20 17:34  elliptic12i.m
     文件        2971  2005-10-08 14:29  elliptic3.m
     文件        2537  2009-06-20 17:35  inversenomeq.m
     文件        4837  2009-06-20 17:35  jacobiThetaEta.m
     文件         968  2009-06-20 17:36  nomeq.m
     文件        2569  2009-06-20 17:36  theta.m
     文件        1312  2014-02-12 12:51  license.txt