www.gusucode.com > elfun18工具箱matlab源码程序 > elfun18/elfun18v1_3/examples/test_verification/test_BF_ident_jacobi.m
% Identities from % Byrd, P.F., Friedman, M.D., 1971. "Handbook of elliptic integrals for % engineers and scientists", 2d ed. Springer-Verlag, Berlin, New York,. cfun = {... 'jcd(x,k) - jcn(x,k)/jdn(x,k)', ... 'jcs(x,k) - jcn(x,k)/jsn(x,k)', ... 'jdc(x,k) - jdn(x,k)/jcn(x,k)', ... 'jds(x,k) - jdn(x,k)/jsn(x,k)', ... 'jsc(x,k) - jsn(x,k)/jcn(x,k)', ... 'jsd(x,k) - jsn(x,k)/jdn(x,k)', ... 'jnc(x,k) - 1/jcn(x,k)', ... 'jnd(x,k) - 1/jdn(x,k)', ... 'jns(x,k) - 1/jsn(x,k)', ... 'jsn(x,k)^2 + jcn(x,k)^2 - 1',... % 121.00 'jdn(x,k)^2 + k^2*jsn(x,k)^2 - 1',... % 121.00 'jdn(x,k)^2 - k^2*jcn(x,k)^2 - kc^2',... % 121.00 'jcn(x,k)^2 + k^2*jsn(x,k)^2 - jdn(x,k)^2',... % 121.00 'jam(-x,k) + jam(x,k)',... % 122.00 'jsn(-x,k) + jsn(x,k)',... 'jcn(-x,k) - jcn(x,k)',... 'jdn(-x,k) - jdn(x,k)',... 'jsc(-x,k) + jsc(x,k)',... 'jam(0,k)',... % 122.01 'jsn(0,k)',... 'jcn(0,k) - 1',... 'jdn(0,k) - 1',... 'jsc(0,k)',... 'jam(elK(k),k) - pi/2',... % 122.02 'jsn(elK(k),k) - 1',... 'jcn(elK(k),k)' ,... 'jdn(elK(k),k) - kc',... 'jam(x + elK(k),k) - atan(kc*jsc(x,k)) - pi/2',... % 122.03 'jsn(x + elK(k),k) - jcd(x,k)',... 'jcn(x + elK(k),k) + kc*jsd(x,k)',... 'jdn(x + elK(k),k) - kc*jnd(x,k)',... 'jsc(x + elK(k),k) + jcs(x,k)/kc',... 'jam(x + 2*elK(k),k) - jam(x,k) - pi',... % 122.04 'jsn(x + 2*elK(k),k) + jsn(x,k)',... % 122.04 'jcn(x + 2*elK(k),k) + jcn(x,k)',... % 122.04 'jdn(x + 2*elK(k),k) - jdn(x,k)',... % 122.04 'jsc(x + 2*elK(k),k) - jsc(x,k)',... % 122.04 'jam(x + 3*elK(k),k) - 3*pi/2 - atan(kc*jsc(x,k))',... % 122.05 'jsn(x + 3*elK(k),k) + jcd(x,k)',... % 122.05 'jcn(x + 3*elK(k),k) - kc*jsd(x,k)',... % 122.05 'jdn(x + 3*elK(k),k) - kc*jnd(x,k)',... % 122.05 'jsc(x + 3*elK(k),k) + jcs(x,k)/kc',... % 122.05 'jam(x + 4*elK(k),k) - jam(x,k) - 2*pi', ... % 122.06 'jsn(x + 4*elK(k),k) - jsn(x,k)', ... % 122.06 'jcn(x + 4*elK(k),k) - jcn(x,k)', ... % 122.06 'jdn(x + 4*elK(k),k) - jdn(x,k)', ... % 122.06 'jsc(x + 4*elK(k),k) - jsc(x,k)',... % 122.06 'jsn(x,0) - sin(x)',... % 122.08 'jcn(x,0) - cos(x)',... 'jdn(x,0) - 1',... 'jsc(x,0) - tan(x)',... 'jsn(x,1) - tanh(x)',... % 122.09 'jcn(x,1) - sech(x)',... 'jdn(x,1) - sech(x)',... 'jsc(x,1) - sinh(x)',... 'jsn(elK(k)/2,k) - 1/sqrt(1 + kc)',... % 122.10 'jcn(elK(k)/2,k) - sqrt(kc/(1 + kc))',... % 122.10 'jdn(elK(k)/2,k) - sqrt(kc)',... % 122.10 'jsc(elK(k)/2,k) - 1/sqrt(kc)',... % 122.10 'jsn(3*elK(k)/2,k) - 1/sqrt(1 + kc)',... % 122.11 'jcn(3*elK(k)/2,k) + sqrt(kc/(1 + kc))',... % 122.11 'jdn(3*elK(k)/2,k) - sqrt(kc)',... % 122.11 'jsc(3*elK(k)/2,k) + 1/sqrt(kc)',... % 122.10 '(1 - jcn(2*x,k))/(1 + jcn(2*x,k)) - jsc(x,k)^2*jdn(x,k)^2', ... %129.04 '(1 - jdn(2*x,k))/(1 + jdn(2*x,k)) - k^2*jsn(x,k)^2*jcd(x,k)^2' ... %129.04 }; bplot = false; nrun = 10000; ntest = length(cfun); NMX = 10; fmt = strcat(repmat('%10.2g',1,8),'\n'); fmt1 = strcat(repmat('%10s',1,8),'\n'); clin = repmat('-',1,90); fprintf(' Consistency tests for elliptic integrals functions \n') fprintf(' Identities from Byrd-Friedman -- Handbook of Elliptic Integrals\n') fprintf(' Modulus -1<=k<=1 \n') fprintf('\n'); fprintf('number of runs %g\n',nrun); fprintf('fail = aerr>sqrt(eps)\n'); for nt = 1:ntest kk = zeros(nrun,1); zz = zeros(nrun,1); nfail = 0; f = zeros(nrun,1); fun = str2func(strcat('@(x,k,kc)',cfun{nt})); fprintf('\ntest %d %s = 0\n',nt, cfun{nt}) fprintf('%s\n',clin); fprintf(fmt1,'min |x|','max |x|','min |k|','max |k|',... 'MAE/eps','RME/eps','NDIG','nfail%'); fprintf('%s\n',clin); rng('shuffle'); for nn = 1:NMX N = 2^nn; for n = 1:nrun k = randi([-N, N]); z = randi([-N, N]); k = (-1)^(n+1)*rand()*2^k; z = (-1)^n*rand()*2^z; if abs(k) > 1 k = 1/k; end kc = sqrt((1 - k)*(1 + k)); f(n) = fun(z,k,kc); if isnan(f(n)) || isinf(f(n)) f(n) = 0; continue end f(n) = abs(f(n)); if f(n) > sqrt(eps) nfail = nfail + 1; end kk(n) = k; zz(n) = z; end fprintf(fmt,... min(abs(zz)),max(abs(zz)),min(abs(kk)), max(abs(kk)),... max(f)/eps,rms(f)/eps,log10(max(f)),... nfail/nrun*100); if nfail > 10 break end end fprintf('%70s\n',repmat('-',1,90)); if ~bplot continue end figure(nt) clf ctitle = strcat('Test ',num2str(nt)); mstp = 100; subplot(1,2,1) hold on title(ctitle) scatter(log10(kk(1:mstp:nrun)),log10(f(1:mstp:nrun))) xlabel('log10(k)') ylabel('log10(aerr)') grid on hold off subplot(1,2,2) hold on scatter(log10(zz(1:mstp:nrun)),log10(f(1:mstp:nrun))) xlabel('log10(u)') ylabel('log10(aerr)') grid on hold off end