標籤（空格分隔）： matlab 積分 數值 計算

matlab數值積分

1 Gauss-Hermite積分

1.1 測試Gauss-Hermite積分函式gaussHermiteIntegral()

clc,clear
format long
[email protected](x)(exp(-x.^4))
% % % % [email protected](x)((1/sqrt(2*pi))*exp(-x.^2/2)) %標準正太分佈
% % area=[];nn=[];
% % for n=1:10
% %     nn=[nn,n];
% %     area=[area,gaussHermiteIntegral(f,-inf,inf,n)];
% % end,nn,area
% % darea=abs(diff(area));
% % %% 系統計算的積分精確值
% % 
 
[email protected](x)(f(-x))
% % area_system=quadgk(f,0,inf,'RelTol',1e-8,'AbsTol',1e-12)+...
% %             quadgk(f_,0,inf,'RelTol',1e-8,'AbsTol',1e-12)
% % d_my_system=abs(area-area_system);
% % log_darea=round(log10(darea)*10)
% % log_d_my_system=round(log10(d_my_system)*10)
%    1.812804954110955
n=5
[area,gErrer]=gaussHermiteIntegral(f,-inf,inf,n,true)
%% 從以下資料看（前向差商，系統誤差）：前向差商公式完全能夠作為事後誤差估計
% log_darea =   -19   -28   -36   -45   -52   -60   -69   -90   -81
% log_d_my_system = -19   -29   -37   -45   -53   -61   -69   -82   -83   -86
%% 從以下資料看：Gauss Hermite積分的精確度e隨積分節點個數n呈現關係：log10(e)=n/10
        %因此採用10n個節點的gauss積分作為程式中的誤差階數
% f =  @(x)(exp(-x.^4))
% n =100
% area = 1.812804956399165
% area_system = 1.812804954110955
% d_my_system = 2.288210509959754e-09
% n=41
% d_my_system = 2.403635556458283e-05
% n=20
% d_my_system =   0.001299981246872
% n=10
% d_my_system =  0.014102477783775
 

執行結果：

f = 
    @(x)(exp(-x.^4))
n =
     5
area =
   1.812809811622279
gErrer =
     3.459687060525241e-05

1.2 gaussHermiteIntegral()

function [area,gErrer]=gaussHermiteIntegral(func,a,b,n,withError)
%高斯積分
%理論依據：數值分析方法 奚梅成
%func 被積函式；
%[a,b] 積分割槽間
%n 高斯積分的階數
%%
%name：鄧能財     Date: 2013/12/23
%% 預設引數
if nargin<5 withError=false; gErrer=0; end
% 查錯
assert(a==-inf && b==inf,'積分割槽間必須為(-inf,inf)!');
assert(ismember(n,[1:10]),'積分精確度數量級必須在1~10之間!');
%% 函式變換
 
[email protected](x)(exp(x^2)*func(x));
%% 資料：10,20,30,...,maxn*10階高斯積分的積分系數
maxn=10; %   (±積分節點node，積分系數coef)，‘±’被省略
data= struct('node',cell(1,maxn),'coef',cell(1,maxn));
%data()= struct('node',{[]},'coef',{[]});
data(1)= struct('node',{[3.4361591188377e+00   
2.5327316742328e+00   
1.7566836492999e+00   
1.0366108297895e+00   
3.4290132722370e-01   ]},'coef',{[7.6404328552326e-06   
1.3436457467812e-03   
3.3874394455481e-02   
2.4013861108231e-01   
6.1086263373533e-01   ]});
data(2)= struct('node',{[5.3874808900112e+00   
4.6036824495507e+00   
3.9447640401156e+00   
3.3478545673832e+00   
2.7888060584281e+00   
2.2549740020893e+00   
1.7385377121166e+00   
1.2340762153953e+00   
7.3747372854539e-01   
2.4534070830090e-01   ]},'coef',{[2.2293936455341e-13   
4.3993409922731e-10   
1.0860693707693e-07   
7.8025564785321e-06   
2.2833863601635e-04   
3.2437733422379e-03   
2.4810520887464e-02   
1.0901720602002e-01   
2.8667550536283e-01   
4.6224366960061e-01   ]});
data(3)= struct('node',{[6.8633452935299e+00   
6.1382792201239e+00   
5.5331471515675e+00   
4.9889189685899e+00   
4.4830553570925e+00   
4.0039086038612e+00   
3.5444438731553e+00   
3.0999705295864e+00   
2.6671321245356e+00   
2.2433914677615e+00   
1.8267411436037e+00   
1.4155278001982e+00   
1.0083382710467e+00   
6.0392105862555e-01   
2.0112857654887e-01   ]},'coef',{[2.9082547163104e-21   
2.8103336037565e-17   
2.8786070806795e-14   
8.1061862974598e-12   
9.1785804243669e-10   
5.1085224507758e-08   
1.5790948873248e-06   
2.9387252289230e-05   
3.4831012431868e-04   
2.7379224730677e-03   
1.4703829704827e-02   
5.5144176870234e-02   
1.4673584754089e-01   
2.8013093083921e-01   
3.8639488954181e-01   ]});
data(4)= struct('node',{[8.0987611392509e+00   
7.4115825314855e+00   
6.8402373052494e+00   
6.3282553512201e+00   
5.8540950560304e+00   
5.4066542479701e+00   
4.9792609785453e+00   
4.5675020728444e+00   
4.1682570668325e+00   
3.7792067534352e+00   
3.3985582658596e+00   
3.0248798839013e+00   
2.6569959984429e+00   
2.2939171418751e+00   
1.9347914722823e+00   
1.5788698949316e+00   
1.2254801090463e+00   
8.7400661235709e-01   
5.2387471383228e-01   
1.7453721459758e-01   ]},'coef',{[2.5910578133142e-29   
8.5440614848689e-25   
2.5675935123395e-21   
1.9891810183765e-18   
6.0083587891597e-16   
8.8057076445574e-14   
7.1565280526724e-12   
3.5256207913729e-10   
1.1212360832276e-08   
2.4111441636703e-07   
3.6315761506931e-06   
3.9369339810925e-05   
3.1385359454133e-04   
1.8714968295980e-03   
8.4608880082581e-03   
2.9312565536172e-02   
7.8474605865404e-02   
1.6337873271327e-01   
2.6572825187738e-01   
3.3864327742559e-01   ]});
data(5)= struct('node',{[9.1824069581293e+00   
8.5227710309178e+00   
7.9756223682056e+00   
7.4864094298642e+00   
7.0343235097706e+00   
6.6086479738554e+00   
6.2029525192747e+00   
5.8129946754204e+00   
5.4357860872249e+00   
5.0691175849172e+00   
4.7112936661690e+00   
4.3609731604546e+00   
4.0170681728581e+00   
3.6786770625153e+00   
3.3450383139379e+00   
3.0154977695745e+00   
2.6894847022677e+00   
2.3664939042987e+00   
2.0460719686864e+00   
1.7278065475159e+00   
1.4113177548983e+00   
1.0962511289577e+00   
7.8227172955461e-01   
4.6905905667824e-01   
1.5630254688947e-01   ]},'coef',{[1.8283576319454e-37   
1.6732639420678e-32   
1.2152047767042e-28   
2.1376568843541e-25   
1.4170933459735e-22   
4.4709843726175e-20   
7.7423829723092e-18   
8.0942618938896e-16   
5.4659440315878e-14   
2.5066555238998e-12   
8.1118773649361e-11   
1.9090405438119e-09   
3.3467934040213e-08   
4.4570299668179e-07   
4.5816827079555e-06   
3.6840190537807e-05   
2.3426989210925e-04   
1.1890117817496e-03   
4.8532638261719e-03   
1.6031941068412e-02   
4.3079159156766e-02   
9.4548935477086e-02   
1.7003245567716e-01   
2.5113085633200e-01   
3.0508512920440e-01   ]});
data(6)= struct('node',{[1.0159109246180e+01   
9.5209036770133e+00   
8.9923980014049e+00   
8.5205692841176e+00   
8.0851886542490e+00   
7.6758399375049e+00   
7.2862765943956e+00   
6.9123815321893e+00   
6.5512591670629e+00   
6.2007735579934e+00   
5.8592901963942e+00   
5.5255210861387e+00   
5.1984265345763e+00   
4.8771500774732e+00   
4.5609737579358e+00   
4.2492864359560e+00   
3.9415607339262e+00   
3.6373358761707e+00   
3.3362046535476e+00   
3.0378033382307e+00   
2.7418037480697e+00   
2.4479069023077e+00   
2.1558378712292e+00   
1.8653415312330e+00   
1.5761790119750e+00   
1.2881246748689e+00   
1.0009634995607e+00   
7.1448878167258e-01   
4.2850006422063e-01   
1.4280123870344e-01   ]},'coef',{[3.7989555024173e-45   
3.2023914909894e-40   
3.9179142693654e-36   
1.3335820885372e-32   
1.7163780212922e-29   
1.0294248187970e-26   
3.3457568162026e-24   
6.5125663921413e-22   
8.1536403865806e-20   
6.9232479198565e-18   
4.1524441111527e-16   
1.8166245761912e-14   
5.9484305160636e-13   
1.4889573490645e-11   
2.8993590128044e-10   
4.4568227752248e-09   
5.4755546192770e-08   
5.4335161342049e-07   
4.3942869362670e-06   
2.9187419041555e-05   
1.6027733468185e-04   
7.3177355696551e-04   
2.7913248289531e-03   
8.9321783603078e-03   
2.4061272766109e-02   
5.4718970932183e-02   
1.0529876369779e-01   
1.7177615691888e-01   
2.3786890495866e-01   
2.7985311752283e-01   ]});
data(7)= struct('node',{[1.1055240743138e+01   
1.0434459776321e+01   
9.9210064725726e+00   
9.4631239608462e+00   
9.0410623133341e+00   
8.6446518655005e+00   
8.2677951152318e+00   
7.9064749419688e+00   
7.5578680649896e+00   
7.2198940023019e+00   
6.8909628903904e+00   
6.5698241787567e+00   
6.2554707188783e+00   
5.9470752244714e+00   
5.6439466223452e+00   
5.3454991415733e+00   
5.0512298489716e+00   
4.7607019531714e+00   
4.4735321501382e+00   
4.1893808634926e+00   
3.9079445988681e+00   
3.6289498686373e+00   
3.3521483007724e+00   
3.0773126524498e+00   
2.8042335229303e+00   
2.5327166122907e+00   
2.2625804097899e+00   
1.9936542226231e+00   
1.7257764756010e+00   
1.4587932269484e+00   
1.1925568563654e+00   
9.2692488970662e-01   
6.6175893081257e-01   
3.9692367564299e-01   
1.3228598727032e-01   ]},'coef',{[5.9928891485445e-50   
1.3088219454484e-44   
1.9775381832792e-41   
5.1459796204458e-40   
1.2902253872744e-36   
1.3859983185653e-33   
7.5869522154891e-31   
2.5172152573746e-28   
5.3429813146915e-26   
7.6615884430059e-24   
7.7454741231880e-22   
5.7095124891098e-20   
3.1533210598112e-18   
1.3342434683070e-16   
4.4061175397096e-15   
1.1534952560119e-13   
2.4259623022220e-12   
4.1457910332392e-11   
5.8137344484274e-10   
6.7473312881327e-09   
6.5292934189487e-08   
5.3024856335052e-07   
3.6344904047183e-06   
2.1131233267306e-05   
1.0466994285007e-04   
4.4339999663333e-04   
1.6117353953559e-03   
5.0416367476379e-03   
1.3605327560619e-02   
3.1741278190051e-02   
6.4133640985369e-02   
1.1238816501982e-01   
1.7101015293604e-01   
2.2612844383734e-01   
2.5999310620316e-01   ]});
data(8)= struct('node',{[1.1887863560471e+01   
1.1281694164989e+01   
1.0780796472469e+01   
1.0334491039538e+01   
9.9234351145692e+00   
9.5376679163021e+00   
9.1712174787265e+00   
8.8201501286611e+00   
8.4817021873216e+00   
8.1538381574635e+00   
7.8350036783293e+00   
7.5239772908069e+00   
7.2197764697589e+00   
6.9215953729258e+00   
6.6287620808576e+00   
6.3407083213300e+00   
6.0569474734510e+00   
5.7770582280616e+00   
5.5006722123152e+00   
5.2274644551007e+00   
4.9571459285134e+00   
4.6894576329528e+00   
4.4241658477652e+00   
4.1610582741258e+00   
3.8999408693874e+00   
3.6406352232272e+00   
3.3829763625124e+00   
3.1268108983731e+00   
2.8719954485287e+00   
2.6183952824718e+00   
2.3658831480742e+00   
2.1143382465070e+00   
1.8636453287449e+00   
1.6136938918511e+00   
1.3643774570540e+00   
1.1155929146033e+00   
8.6723992270397e-01   
6.1922034962757e-01   
3.7143774948304e-01   
1.2379686317313e-01   ]},'coef',{[3.8685383727931e-49   
2.0252403350315e-45   
5.6114941460751e-43   
2.7312591323502e-40   
8.0059947621974e-39   
4.8548677202834e-38   
1.6512006592794e-36   
6.3922427517268e-35   
1.8964852968049e-32   
4.3114531467056e-30   
6.8829012692744e-28   
7.9853336688990e-26   
6.9341381005532e-24   
4.6134698831073e-22   
2.3978402197188e-20   
9.8965290426435e-19   
3.2891620153310e-17   
8.9094610374559e-16   
1.9875208348038e-14   
3.6848498201918e-13   
5.7232797184518e-12   
7.4997199146582e-11   
8.3430559260916e-10   
7.9229232891485e-09   
6.4544629501197e-08   
4.5305355658615e-07   
2.7507246274012e-06   
1.4496444191130e-05   
6.6517412418152e-05   
2.6647788669334e-04   
9.3431768696117e-04   
2.8731975330897e-03   
7.7640394053997e-03   
1.8465751154428e-02   
3.8708599914838e-02   
7.1600747691843e-02   
1.1698073380430e-01   
1.6893796927590e-01   
2.1577494627433e-01   
2.4383585380721e-01   ]});
data(9)= struct('node',{[1.2668764375434e+01   
1.2075130816911e+01   
1.1584960252805e+01   
1.1148508777694e+01   
1.0746787108514e+01   
1.0370015493271e+01   
1.0012330403392e+01   
9.6698698743375e+00   
9.3399210724884e+00   
9.0204865095226e+00   
8.7100414635652e+00   
8.4073884130035e+00   
8.1115647551405e+00   
7.8217816701739e+00   
7.5373821269924e+00   
7.2578111510466e+00   
6.9825942253569e+00   
6.7113212484547e+00   
6.4436343875168e+00   
6.1792187234935e+00   
5.9177949371442e+00   
5.6591135131207e+00   
5.4029500908349e+00   
5.1491016937770e+00   
4.8973836402140e+00   
4.6476269884214e+00   
4.3996764055654e+00   
4.1533883754886e+00   
3.9086296798929e+00   
3.6652761017604e+00   
3.4232113106648e+00   
3.1823258978579e+00   
2.9425165353484e+00   
2.7036852380959e+00   
2.4657387122777e+00   
2.2285877756048e+00   
1.9921468380504e+00   
1.7563334332534e+00   
1.5210677923813e+00   
1.2862724534456e+00   
1.0518719000410e+00   
8.1779222425490e-01   
5.8396080911088e-01   
3.5030602639582e-01   
1.1675694609164e-01   ]},'coef',{[1.3926602246895e-51   
1.9568360107978e-47   
1.0639144183847e-44   
7.9089094484501e-43   
7.1933465481042e-41   
3.4996773126588e-39   
9.2942987391649e-38   
5.0159398402133e-38   
3.1216254599340e-37   
2.8172065374851e-36   
3.9144464842405e-34   
5.9935373474292e-32   
7.7735858553266e-30   
7.7182478105170e-28   
5.9828467035793e-26   
3.6826786261206e-24   
1.8269175754861e-22   
7.3963427902401e-21   
2.4705506323895e-19   
6.8737797964624e-18   
1.6064742921236e-16   
3.1773129595293e-15   
5.3533724860157e-14   
7.7293031086162e-13   
9.6137828692418e-12   
1.0350233726258e-10   
9.6863849248741e-10   
7.9104166149095e-09   
5.6567852900933e-08   
3.5533045829498e-07   
1.9661408494659e-06   
9.6077760741792e-06   
4.1557921102543e-05   
1.5944214368481e-04   
5.4359224952963e-04   
1.6496055361812e-03   
4.4622998595035e-03   
1.0773786624642e-02   
2.3243228173755e-02   
4.4849988132390e-02   
7.7467627605513e-02   
1.1985655590239e-01   
1.6619497260996e-01   
2.0661430370392e-01   
2.3035797018195e-01   ]});
data(10)= struct('node',{[1.3406487338145e+01   
1.2823799749488e+01   
1.2342964222860e+01   
1.1915061943114e+01   
1.1521415400787e+01   
1.1152404385585e+01   
1.0802260753685e+01   
1.0467185421343e+01   
1.0144509941293e+01   
9.8322698077780e+00   
9.5289658233901e+00   
9.2334208902192e+00   
8.9446892173255e+00   
8.6619961681345e+00   
8.3846969404163e+00   
8.1122473111628e+00   
7.8441823844608e+00   
7.5801008078575e+00   
7.3196528223045e+00   
7.0625310602489e+00   
6.8084633528588e+00   
6.5572070319215e+00   
6.3085443611121e+00   
6.0622788326143e+00   
5.8182321352035e+00   
5.5762416493299e+00   
5.3361583601384e+00   
5.0978451050891e+00   
4.8611750917912e+00   
4.6260306357872e+00   
4.3923020786827e+00   
4.1598868551310e+00   
3.9286886834277e+00   
3.6986168593185e+00   
3.4695856364186e+00   
3.2415136796310e+00   
3.0143235803312e+00   
2.7879414239820e+00   
2.5622964023726e+00   
2.3373204639069e+00   
2.1129479963712e+00   
1.8891155374270e+00   
1.6657615087415e+00   
1.4428259702159e+00   
1.2202503912190e+00   
9.9797743609810e-01   
7.7595076154015e-01   
5.5411482359162e-01   
3.3241469234223e-01   
1.1079587242244e-01   ]},'coef',{[4.1890734835153e-53   
6.6685952840498e-49   
1.7149047609884e-45   
4.6329593470698e-43   
2.3803897130070e-41   
3.1234305446425e-39   
2.4962135456919e-38   
2.5988685504385e-37   
1.3109926053898e-36   
1.4501664855704e-36   
5.1295692004930e-37   
1.3477869587572e-35   
2.8227154103471e-35   
5.9582141086172e-34   
7.7461504141821e-32   
7.1102466014413e-30   
5.0387971358429e-28   
2.9171459549279e-26   
1.3948232985713e-24   
5.5610285249618e-23   
1.8649975920043e-21   
5.3023160356162e-20   
1.2868329542238e-18   
2.6824921883091e-17   
4.8298353231531e-16   
7.5488968761282e-15   
1.0288749372616e-13   
1.2278785143771e-12   
1.2879038257424e-11   
1.1913006349369e-10   
9.7479212538641e-10   
7.0758572838838e-09   
4.5681275084851e-08   
2.6290974837540e-07   
1.3517971591104e-06   
6.2215248177779e-06   
2.5676159384549e-05   
9.5171627785510e-05   
3.1729197104330e-04   
9.5269218854862e-04   
2.5792732600591e-03   
6.3030002856081e-03   
1.3915665220232e-02   
2.7779127385933e-02   
5.0175812677429e-02   
8.2051827391224e-02   
1.2153798684410e-01   
1.6313003050278e-01   
1.9846285025419e-01   
2.1889262958744e-01   ]});
%% 計算
area=0;
noden=data(n).node;
coefn=data(n).coef;
k=n;
n=10*n;  %n轉換為表示Gauss積分的積分節點或積分系數的個數
if n/2-floor(n/2)==1/2 %當階數為奇數
    for i=1:floor(n/2)
        area=area+coefn(i)*(expfunc(noden(i))+expfunc(-noden(i)));
    end
    i=floor(n/2)+1;
    area=area+coefn(i)*expfunc(noden(i));
elseif n/2-floor(n/2)==0 %當階數為偶數
    for i=1:n/2
        area=area+coefn(i)*(expfunc(noden(i))+expfunc(-noden(i)));
    end
end
%% 誤差:前向（n-1階）差商作為事後誤差估計
if withError
    area_=gaussHermiteIntegral(func,a,b,k-1,false);
    gErrer=abs(area-area_);
end
end

2 Gauss-Laguerre積分

2.1 測試Gauss-Laguerre積分函式gaussLaguerreIntegral()

clc,clear
format long
% [email protected](x)((x+1).^(-2))
% [email protected](x)(exp(-x.^2))
% [email protected](x)(exp(-x))
% [email protected](x)(exp(-x.^4))
[email protected](x)((1/sqrt(2*pi))*exp(-x.^2/2)) %標準正太分佈
area=[];nn=[];
for n=1:10
    nn=[nn,n];
    area=[area,gaussLaguerreIntegral(f,0,inf,n)];
end,nn
area_str=printArray(area)
darea=abs(diff(area));
%% 系統計算的積分精確值
% area_system=quadgk(f,0,inf,'RelTol',1e-8,'AbsTol',1e-12)
d_my_system=abs(area-.5);
d_my_system_str=printArray(d_my_system)
log_darea=round(log10(darea)*10);
log_d_my_system=round(log10(d_my_system)*10);
log_darea_str=printArray(log_darea)
log_d_my_system_str=printArray(log_d_my_system)
%% 測試誤差項
n=5
[area,gErrer]=gaussLaguerreIntegral(f,0,inf,n,true)

執行結果：

f = 
    @(x)((1/sqrt(2*pi))*exp(-x.^2/2))
n =
     5
area =
   0.500000012302732
gErrer =
     9.329705445981773e-07

2.2 gaussLaguerreIntegral()

function [area,gErrer]=gaussLaguerreIntegral(func,a,b,n,withError)
%高斯積分
%理論依據：數值分析方法 奚梅成
%func 被積函式；
%[a,b] 積分割槽間
%n 高斯積分的階數
%%
%name：鄧能財     Date: 2013/12/21
%% 預設引數
if nargin<5 withError=false; gErrer=0; end
%% 查錯
maxn=10;
% assert(a==0 && b==inf,'積分割槽間必須為[0,inf)!');
% assert(ismember(n,[1:maxn]),['積分精確度數量級必須在1~',num2str(maxn),'之間!']);
%% 函式變換
[email protected](x)(exp(x)*func(x));
%% 資料：h,2h,3h,...,maxn*h階高斯積分的積分系數
h=5;
    %   (積分節點node，積分系數coef)
data= struct('node',cell(1,maxn),'coef',cell(1,maxn));
%data()= struct('node',{[]},'coef',{[]});
data(1)= struct('node',{[2.6356031971814e-01
    1.4134030591065e+00
    3.5964257710407e+00
    7.0858100058588e+00
    1.2640800844276e+01   ]},'coef',{[5.2175561058281e-01
    3.9866681108317e-01
    7.5942449681707e-02
    3.6117586799220e-03
    2.3369972385776e-05   ]});
data(2)= struct('node',{[1.3779347054049e-01
    7.2945454950317e-01
    1.8083429017403e+00
    3.4014336978549e+00
    5.5524961400638e+00
    8.3301527467645e+00
    1.1843785837900e+01
    1.6279257831378e+01
    2.1996585811981e+01
    2.9920697012274e+01   ]},'coef',{[3.0844111576502e-01
    4.0111992915527e-01
    2.1806828761181e-01
    6.2087456098678e-02
    9.5015169751811e-03
    7.5300838858754e-04
    2.8259233495996e-05
    4.2493139849627e-07
    1.8395648239796e-09
    9.9118272196090e-13   ]});
data(3)= struct('node',{[9.3307812017283e-02
    4.9269174030189e-01
    1.2155954120710e+00
    2.2699495262037e+00
    3.6676227217514e+00
    5.4253366274136e+00
    7.5659162266131e+00
    1.0120228568019e+01
    1.3130282482176e+01
    1.6654407708330e+01
    2.0776478899449e+01
    2.5623894226729e+01
    3.1407519169754e+01
    3.8530683306486e+01
    4.8026085572686e+01   ]},'coef',{[2.1823488594009e-01
    3.4221017792288e-01
    2.6302757794168e-01
    1.2642581810593e-01
    4.0206864921001e-02
    8.5638778036118e-03
    1.2124361472143e-03
    1.1167439234425e-04
    6.4599267620229e-06
    2.2263169070963e-07
    4.2274303849794e-09
    3.9218972670411e-11
    1.4565152640731e-13
    1.4830270511133e-16
    1.6005949062111e-20   ]});
data(4)= struct('node',{[7.0539889691990e-02
    3.7212681800161e-01
    9.1658210248327e-01
    1.7073065310283e+00
    2.7491992553094e+00
    4.0489253138509e+00
    5.6151749708616e+00
    7.4590174536711e+00
    9.5943928695811e+00
    1.2038802546964e+01
    1.4814293442631e+01
    1.7948895520519e+01
    2.1478788240285e+01
    2.5451702793187e+01
    2.9932554631701e+01
    3.5013434240479e+01
    4.0833057056729e+01
    4.7619994047347e+01
    5.5810795750064e+01
    6.6524416525616e+01   ]},'coef',{[1.6874680185112e-01
    2.9125436200607e-01
    2.6668610286700e-01
    1.6600245326951e-01
    7.4826064668792e-02
    2.4964417309283e-02
    6.2025508445722e-03
    1.1449623864769e-03
    1.5574177302781e-04
    1.5401440865225e-05
    1.0864863665180e-06
    5.3301209095567e-08
    1.7579811790506e-09
    3.7255024025123e-11
    4.7675292515782e-13
    3.3728442433624e-15
    1.1550143395004e-17
    1.5395221405823e-20
    5.2864427255692e-24
    1.6564566124990e-28   ]});
data(5)= struct('node',{[5.6704775452705e-02
    2.9901089858699e-01
    7.3590955543502e-01
    1.3691831160352e+00
    2.2013260537215e+00
    3.2356758035580e+00
    4.4764966150738e+00
    5.9290837627004e+00
    7.5998993099567e+00
    9.4967492209324e+00
    1.1629014911779e+01
    1.4007957976545e+01
    1.6647125597289e+01
    1.9562898011469e+01
    2.2775241986835e+01
    2.6308772390969e+01
    3.0194291163316e+01
    3.4471097571922e+01
    3.9190608803937e+01
    4.4422349336162e+01
    5.0264574993834e+01
    5.6864967173940e+01
    6.4466670615954e+01
    7.3534234792100e+01
    8.5260155562496e+01   ]},'coef',{[1.3752601422934e-01
    2.5164527376491e-01
    2.5617600280975e-01
    1.8621549036244e-01
    1.0319984810752e-01
    4.4714161129934e-02
    1.5305232886396e-02
    4.1524146328771e-03
    8.9209907325968e-04
    1.5115601916424e-04
    2.0065531801933e-05
    2.0677743964319e-06
    1.6346520222911e-07
    9.7660150621244e-09
    4.3277207941849e-10
    1.3896009633895e-11
    3.1389227925400e-13
    4.8026148226043e-15
    4.7358853648073e-17
    2.8142053798431e-19
    9.1649543959912e-22
    1.4189400094973e-24
    8.2736519440991e-28
    1.1688817115428e-31
    1.3158315000591e-36   ]});data(6)= struct('node',{[4.7407180540804e-02   
2.4992391675316e-01   
6.1483345439277e-01   
1.1431958256661e+00   
1.8364545546226e+00   
2.6965218745572e+00   
3.7258145077795e+00   
4.9272937658499e+00   
6.3045155909651e+00   
7.8616932933703e+00   
9.6037759854793e+00   
1.1536546597956e+01   
1.3666744693064e+01   
1.6002221188981e+01   
1.8552134840143e+01   
2.1327204321783e+01   
2.4340035764533e+01   
2.7605554796781e+01   
3.1141586701111e+01   
3.4969652008249e+01   
3.9116084949068e+01   
4.3613652908485e+01   
4.8503986163804e+01   
5.3841385406508e+01   
5.9699121859235e+01   
6.6180617794438e+01   
7.3441238595560e+01   
8.1736810506728e+01   
9.1556466522537e+01   
1.0415752443106e+02   ]},'coef',{[1.1604408602039e-01   
2.2085112475070e-01   
2.4139982758787e-01   
1.9463676844642e-01   
1.2372841596688e-01   
6.3678780368988e-02   
2.6860475273381e-02   
9.3380708816042e-03   
2.6806968913369e-03   
6.3512912194088e-04   
1.2390745990689e-04   
1.9828788438953e-05   
2.5893509291317e-06   
2.7409428405359e-07   
2.3328311650249e-08   
1.5807455747823e-09   
8.4274791230481e-11   
3.4851612347541e-12   
1.0990180599189e-13   
2.5883126679993e-15   
4.4378380184999e-17   
5.3659179405818e-19   
4.3939492584790e-21   
2.3114466007031e-23   
7.2761720424752e-26   
1.2409628392253e-28   
9.9460139172742e-32   
3.0921279934282e-35   
3.1849751877368e-39   
1.7664292620495e-43   ]});
data(7)= struct('node',{[4.0729209061714e-02   
2.1468745273514e-01   
5.2801038431934e-01   
9.8138617345908e-01   
1.5757259475775e+00   
2.3122282965132e+00   
3.1923979394931e+00   
4.2180641250297e+00   
5.3914027346738e+00   
6.7149632799157e+00   
8.1917017588376e+00   
9.8250204970927e+00   
1.1618816388907e+01   
1.3577539357308e+01   
1.5706263395763e+01   
1.8010773287835e+01   
2.0497671107147e+01   
2.3174507997799e+01   
2.6049948710370e+01   
2.9133979209544e+01   
3.2438171837678e+01   
3.5976028769894e+01   
3.9763434104800e+01   
4.3819260118607e+01   
4.8166197974019e+01   
5.2831925058016e+01   
5.7850795022134e+01   
6.3266373675385e+01   
6.9135413883648e+01   
7.5534435134730e+01   
8.2571405523583e+01   
9.0408523864401e+01   
9.9312973659441e+01   
1.0979599094491e+02   
1.2317325317538e+02   ]},'coef',{[1.0036223741924e-01   
1.9648238336171e-01   
2.2601456693434e-01   
1.9627166789644e-01   
1.3760074201313e-01   
8.0030400464140e-02   
3.9125771428883e-02   
1.6186726526126e-02   
5.6853255717743e-03   
1.6972232991260e-03   
4.3048506814168e-04   
9.2634053063654e-05   
1.6870264869314e-05   
2.5916721474280e-06   
3.3446109572665e-07   
3.6077854386912e-08   
3.2335690322103e-09   
2.3913229815959e-10   
1.4473323990302e-11   
7.1013552776266e-13   
2.7933695098016e-14   
8.6948402151367e-16   
2.1088459299447e-17   
3.9129108212979e-19   
5.4327218369250e-21   
5.4936982290015e-23   
3.9126007648144e-25   
1.8797445108978e-27   
5.7353329163766e-30   
1.0118890306195e-32   
6.4972453567684e-36   
1.0249186153468e-41   
1.0144829299011e-40   
1.8914407768105e-43   
2.6963083454085e-47   ]});
data(8)= struct('node',{[3.5700394308884e-02   
1.8816228315870e-01   
4.6269428131458e-01   
8.5977296397294e-01   
1.3800108205273e+00   
2.0242091359228e+00   
2.7933693535068e+00   
3.6887026779083e+00   
4.7116411465550e+00   
5.8638508783437e+00   
7.1472479081023e+00   
8.5640170175862e+00   
1.0116634048452e+01   
1.1807892294005e+01   
1.3640933712537e+01   
1.5619285893339e+01   
1.7746905950096e+01   
2.0028232834575e+01   
2.2468249983498e+01   
2.5072560772426e+01   
2.7847480009169e+01   
3.0800145739445e+01   
3.3938657084914e+01   
3.7272245880476e+01   
4.0811492823887e+01   
4.4568603175334e+01   
4.8557763533060e+01   
5.2795611187217e+01   
5.7301863323394e+01   
6.2100179072775e+01   
6.7219370927127e+01   
7.2695158847612e+01   
7.8572802911571e+01   
8.4911231135705e+01   
9.1789874671236e+01   
9.9320808717447e+01   
1.0767244063939e+02   
1.1712230951269e+02   
1.2820184198826e+02   
1.4228004446916e+02   ]},'coef',{[8.8412106190338e-02   
1.7681473909573e-01   
2.1136311701596e-01   
1.9408119531860e-01   
1.4643428242412e-01   
9.3326798435770e-02   
5.0932204361044e-02   
2.3976193015685e-02   
9.7746252467145e-03   
3.4579399930185e-03   
1.0622468938969e-03   
2.8327168532433e-04   
6.5509405003248e-05   
1.3116069073268e-05   
2.2684528787793e-06   
3.3796264822014e-07   
4.3228213222805e-08   
4.7284937709896e-09   
4.4031741042515e-10   
3.4724414847684e-11   
2.3053815448611e-12   
1.2797725979698e-13   
5.8941771706254e-15   
2.2322175699210e-16   
6.8803366293934e-18   
1.7056039068355e-19   
3.3537088111640e-21   
5.1461768068696e-23   
6.0453446308727e-25   
5.3139346735449e-27   
3.3657690232258e-29   
1.2178445191183e-31   
4.2890086271870e-35   
2.3067506710852e-34   
2.3302702800589e-35   
7.9126811843002e-37   
1.4714096228146e-38   
6.8575263257425e-41   
1.0375271306839e-43   
1.2955873193971e-47   ]});
data(9)= struct('node',{[3.1776956867006e-02   
1.6747251657590e-01   
4.1176935189346e-01   
7.6501433412688e-01   
1.2276394084198e+00   
1.8002069428945e+00   
2.4834175141916e+00   
3.2781154527750e+00   
4.1852949337855e+00   
5.2061071301655e+00   
6.3418686461714e+00   
7.5940714118218e+00   
8.9643942355385e+00   
1.0454716247742e+01   
1.2067132515935e+01   
1.3803972172016e+01   
1.5667819467697e+01   
1.7661538267980e+01   
1.9788300611284e+01   
2.2051620115564e+01   
2.4455391203194e+01   
2.7003935367666e+01   
2.9702056032584e+01   
3.2555103986001e+01   
3.5569055951601e+01   
3.8750609641046e+01   
4.2107299705812e+01   
4.5647640501793e+01   
4.9381303694811e+01   
5.3319341780002e+01   
5.7474473058814e+01   
6.1861450326356e+01   
6.6497545839445e+01   
7.1403201447971e+01   
7.6602919389220e+01   
8.2126514283411e+01   
8.8010926374694e+01   
9.4302943609710e+01   
1.0106347092232e+02   
1.0837460196798e+02   
1.1635218467520e+02   
1.2517034745017e+02   
1.3511619015017e+02   
1.4673979649750e+02   
1.6145944790909e+02   ]},'coef',{[7.9003863218578e-02   
1.6064969814262e-01   
1.9788792774995e-01   
1.8978490012961e-01   
1.5161633179829e-01   
1.0374664065632e-01   
6.1658521650479e-02   
3.2073009065685e-02   
1.4666691055305e-02   
5.9109404105920e-03   
2.1021480717925e-03   
6.5997635763383e-04   
1.8287173922646e-04   
4.4687548597413e-05   
9.6188526495207e-06   
1.8207784506475e-06   
3.0249820529759e-07   
4.4003996132627e-08   
5.5894745239333e-09   
6.1800526258555e-10   
5.9265751381681e-11   
4.9096834353441e-12   
3.4975836508939e-13   
2.1317103205750e-14   
1.1051884160140e-15   
4.8425951880898e-17   
1.7802221360264e-18   
5.4452758047095e-20   
1.3728113205342e-21   
2.8219449500093e-23   
4.6711191673266e-25   
6.1366267898675e-27   
6.3012553677705e-29   
4.9199469546978e-31   
2.8553923595618e-33   
2.6513859772935e-35   
1.7302431662816e-35   
3.7028624824552e-36   
1.3728993488078e-37   
1.3069964719048e-38   
3.7845759572530e-40   
6.4931338105968e-44   
1.0904747480551e-45   
5.8188542963473e-48   
7.4815330585462e-52   ]});
data(10)= struct('node',{[2.8630518339375e-02   
1.5088293567693e-01   
3.7094878153490e-01   
6.8909069988105e-01   
1.1056250235399e+00   
1.6209617511025e+00   
2.2356103759152e+00   
2.9501833666418e+00   
3.7653997744058e+00   
4.6820893875593e+00   
5.7011975747849e+00   
6.8237909097945e+00   
8.0510636693908e+00   
9.3843453082584e+00   
1.0825109031549e+01   
1.2374981608757e+01   
1.4035754599830e+01   
1.5809397197845e+01   
1.7698070933350e+01   
1.9704146535462e+01   
2.1830223306578e+01   
2.4079151444412e+01   
2.6454057841253e+01   
2.8958376011937e+01   
3.1595880956623e+01   
3.4370729963090e+01   
3.7287510610551e+01   
4.0351297573586e+01   
4.3567720269995e+01   
4.6943043991603e+01   
5.0484267963130e+01   
5.4199244880169e+01   
5.8096828017249e+01   
6.2187054175689e+01   
6.6481373878445e+01   
7.0992944826619e+01   
7.5737011547727e+01   
8.0731404802478e+01   
8.5997211136463e+01   
9.1559690412534e+01   
9.7449565614851e+01   
1.0370489123669e+02   
1.1037385880764e+02   
1.1751919820311e+02   
1.2522547013347e+02   
1.3361202792273e+02   
1.4285832548925e+02   
1.5326037197260e+02   
1.6538564331668e+02   
1.8069834370921e+02   ]},'coef',{[7.1404726135184e-02   
1.4714860696459e-01   
1.8567162757483e-01   
1.8438538252735e-01   
1.5420116860636e-01   
1.1168536990227e-01   
7.1052885490196e-02   
4.0020276911508e-02   
2.0050623080072e-02   
8.9608512036464e-03   
3.5781124153157e-03   
1.2776171567891e-03   
4.0803024498372e-04   
1.1652883223097e-04   
2.9741704936942e-05   
6.7778425265418e-06   
1.3774795031715e-06   
2.4928861817200e-07   
4.0103543504264e-08   
5.7233317481479e-09   
7.2294342491785e-10   
8.0617101422267e-11   
7.9133931001156e-12   
6.8157366176493e-13   
5.1324267151836e-14   
3.3656247656300e-15   
1.9134763294156e-16   
9.3855895940761e-18   
3.9500701020349e-19   
1.4177501120800e-20   
4.3099605395263e-22   
1.1012486685628e-23   
2.3448759094638e-25   
4.1195425783619e-27   
5.8517301056039e-29   
6.7131788910435e-31   
1.2697958123263e-32   
6.9071714060483e-34   
3.9707498882293e-34   
1.2868077986345e-33   
1.1748415289134e-33   
1.8729029481275e-34   
2.9809110912453e-35   
2.1474665129166e-36   
5.6599382245170e-38   
1.2115127694324e-39   
8.4123926853140e-42   
1.1227629730641e-44   
9.4963283290577e-48   
4.9608486792055e-52   ]});
%% 計算
area=0;
noden=data(n).node;
coefn=data(n).coef;
k=h*n;  %k表示Gauss積分的積分節點或積分系數的個數
for i=1:k
    area=area+coefn(i)*expfunc(noden(i));
end
%% 誤差:前向（n-1階）差商作為事後誤差估計
if withError
    area_=gaussLaguerreIntegral(func,a,b,n-1,false);
    gErrer=abs(area-area_);
end
end

輸出函式：

function str=printArray(a)
str=sprintf('%.13e\t',a);

3 Gauss-Legendre積分W(x)=1

3.1 測試Gauss-Legendre積分W(x)=1函式gaussLegendreIntegral()

clc,clear
format long
% [email protected](x)(exp(x))
% a=1,b=10
[email protected](x)(cos(x))
% a=0,b=pi/2
a=-pi/2,b=pi/2
% a=-pi,b=pi
% a=-1,b=1
disp('各個階數(h*)的gaussLegendreIntegral積分')
h=5
max=5
n_=[2:max];area_=[];gErrer_=[];
for n=2:max
%% 適合區間長為數量級10左右的積分
    [area,gErrer]=gaussLegendreIntegral(f,a,b,n,true);
    area_=[area_,area];gErrer_=[gErrer_,gErrer];
end
n_
area__str=printArray(area_)
gErrer_str=printArray(gErrer_)
area_system=quadl(f,a,b)

執行結果：

f = 
    @(x)(cos(x))
a =
  -1.570796326794897
b =
   1.570796326794897
各個階數(h*)的gaussLegendreIntegral積分
h =
     5
max =
     5
n_ =
     2     3     4     5
area__str =
2.0000000000000e+00	2.0000000000000e+00	2.0000000000000e+00	2.0000000000000e+00	
gErrer_str =
1.1028448909656e-07	3.0864200084579e-14	3.0864200084579e-14	2.7089441800854e-14	
area_system =
   1.999999977471133
>> 

3.2 gaussLegendreIntegral()

function [area,gErrer]=gaussLegendreIntegral(func,a,b,n,withError)
%高斯積分
%理論依據：數值分析方法 奚梅成
%func 被積函式；
%[a,b] 積分割槽間
%n 高斯積分的階數
%% 
%name：鄧能財     Date: 2013/12/21
%% 預設引數
if nargin<5 withError=false; gErrer=0; end
%% 資料：1-8階高斯積分的積分系數
h=5;
maxn=10; %   (±積分節點node，積分系數coef)，‘±’被省略
data= struct('node',cell(1,maxn),'coef',cell(1,maxn));
% data(1)= struct('node',{[0           ]},'coef',{[2]});
% data(2)= struct('node',{[0.5773502692]},'coef',{[1]});
% data(3)= struct('node',{[0.7745966692 0           ]},'coef',{[0.5555555556 0.8888888889]});
% data(4)= struct('node',{[0.8611363116 0.3399810436]},'coef',{[0.3478548451 0.6521451549]});
% data(5)= struct('node',{[0.9061798459 0.5384693101 0           ]},'coef',{[0.2369268851 0.4786286705 0.5688888889]});
% data(6)= struct('node',{[0.9324695142 0.6612093865 0.2386191861]},'coef',{[0.1713244924 0.3607615730 0.4679139346]});
% data(7)= struct('node',{[0.9491079123 0.7415311856 0.4058451514 0           ]},'coef',{[0.1294849662 0.2797053915 0.3818300505 0.4179591837]});
% data(8)= struct('node',{[0.9602898565 0.7966664774 0.5255324099 0.1834346425]},'coef',{[0.1012285363 0.2223810345 0.3137066459 0.3626837834]});
data(1)= struct('node',{[9.0617984593866e-01   
5.3846931010568e-01   
-4.7505928254367e-17   ]},'coef',{[2.3692688505619e-01   
4.7862867049937e-01   
5.6888888888889e-01   ]});
data(2)= struct('node',{[9.7390652851717e-01   
8.6506336668898e-01   
6.7940956829902e-01   
4.3339539412925e-01   
1.4887433898163e-01   ]},'coef',{[6.6671344308688e-02   
1.4945134915058e-01   
2.1908636251598e-01   
2.6926671931000e-01   
2.9552422471475e-01   ]});
data(3)= struct('node',{[9.8799251802049e-01   
9.3727339240071e-01   
8.4820658341043e-01   
7.2441773136017e-01   
5.7097217260854e-01   
3.9415134707756e-01   
2.0119409399743e-01   
1.0426645532607e-16   ]},'coef',{[3.0753241996117e-02   
7.0366047488109e-02   
1.0715922046717e-01   
1.3957067792615e-01   
1.6626920581699e-01   
1.8616100001556e-01   
1.9843148532711e-01   
2.0257824192556e-01   ]});
data(4)= struct('node',{[9.9312859918510e-01   
9.6397192727791e-01   
9.1223442825133e-01   
8.3911697182222e-01   
7.4633190646015e-01   
6.3605368072652e-01   
5.1086700195083e-01   
3.7370608871542e-01   
2.2778585114164e-01   
7.6526521133497e-02   ]},'coef',{[1.7614007139152e-02   
4.0601429800386e-02   
6.2672048334110e-02   
8.3276741576705e-02   
1.0193011981724e-01   
1.1819453196152e-01   
1.3168863844918e-01   
1.4209610931838e-01   
1.4917298647260e-01   
1.5275338713073e-01   ]});
data(5)= struct('node',{[9.9555696979050e-01   
9.7666392145952e-01   
9.4297457122897e-01   
8.9499199787828e-01   
8.3344262876083e-01   
7.5925926303736e-01   
6.7356636847347e-01   
5.7766293024122e-01   
4.7300273144572e-01   
3.6117230580939e-01   
2.4386688372099e-01   
1.2286469261071e-01   
2.3477116549018e-17   ]},'coef',{[1.1393798501026e-02   
2.6354986615032e-02   
4.0939156701307e-02   
5.4904695975835e-02   
6.8038333812357e-02   
8.0140700335000e-02   
9.1028261982964e-02   
1.0053594906705e-01   
1.0851962447426e-01   
1.1485825914571e-01   
1.1945576353578e-01   
1.2224244299031e-01   
1.2317605372672e-01   ]});
data(6)= struct('node',{[9.9689348407465e-01   
9.8366812327975e-01   
9.6002186496831e-01   
9.2620004742927e-01   
8.8256053579205e-01   
8.2956576238277e-01   
7.6777743210483e-01   
6.9785049479332e-01   
6.2052618298924e-01   
5.3662414814202e-01   
4.4703376953809e-01   
3.5270472553088e-01   
2.5463692616789e-01   
1.5386991360858e-01   
5.1471842555318e-02   ]},'coef',{[7.9681924961671e-03   
1.8466468311090e-02   
2.8784707883324e-02   
3.8799192569627e-02   
4.8402672830593e-02   
5.7493156217619e-02   
6.5974229882181e-02   
7.3755974737705e-02   
8.0755895229420e-02   
8.6899787201083e-02   
9.2122522237786e-02   
9.6368737174644e-02   
9.9593420586795e-02   
1.0176238974840e-01   
1.0285265289356e-01   ]});
data(7)= struct('node',{[9.9770656909960e-01   
9.8793576444385e-01   
9.7043761603923e-01   
9.4534514820783e-01   
9.1285426135932e-01   
8.7321912502522e-01   
8.2674989909223e-01   
7.7381025228691e-01   
7.1481450155663e-01   
6.5022436466589e-01   
5.8054534474976e-01   
5.0632277324149e-01   
4.2813754151781e-01   
3.4660155443081e-01   
2.6235294120930e-01   
1.7605106116599e-01   
8.8371343275659e-02   
-8.0123445265982e-18   ]},'coef',{[5.8834334204433e-03   
1.3650828348362e-02   
2.1322979911483e-02   
2.8829260108895e-02   
3.6110115863463e-02   
4.3108422326170e-02   
4.9769370401354e-02   
5.6040816212370e-02   
6.1873671966080e-02   
6.7222285269087e-02   
7.2044794772560e-02   
7.6303457155442e-02   
7.9964942242325e-02   
8.3000593728857e-02   
8.5386653392099e-02   
8.7104446997184e-02   
8.8140530430276e-02   
8.8486794907104e-02   ]});
data(8)= struct('node',{[9.9823770971056e-01   
9.9072623869946e-01   
9.7725994998377e-01   
9.5791681921379e-01   
9.3281280827868e-01   
9.0209880696887e-01   
8.6595950321226e-01   
8.2461223083331e-01   
7.7830565142652e-01   
7.2731825518993e-01   
6.7195668461418e-01   
6.1255388966798e-01   
5.4946712509513e-01   
4.8307580168618e-01   
4.1377920437160e-01   
3.4199409082576e-01   
2.6815218500725e-01   
1.9269758070137e-01   
1.1608407067526e-01   
3.8772417506051e-02   ]},'coef',{[4.5212770985336e-03   
1.0498284531152e-02   
1.6421058381909e-02   
2.2245849194167e-02   
2.7937006980023e-02   
3.3460195282548e-02   
3.8782167974472e-02   
4.3870908185673e-02   
4.8695807635072e-02   
5.3227846983937e-02   
5.7439769099392e-02   
6.1306242492929e-02   
6.4804013456601e-02   
6.7912045815234e-02   
7.0611647391287e-02   
7.2886582395805e-02   
7.4723169057968e-02   
7.6110361900626e-02   
7.7039818164248e-02   
7.7505947978425e-02   ]});
data(9)= struct('node',{[9.9860364518194e-01   
9.9264999844720e-01   
9.8196871503454e-01   
9.6660831039689e-01   
9.4664169099563e-01   
9.2216393671900e-01   
8.9329167175324e-01   
8.6016247596066e-01   
8.2293422050209e-01   
7.8178431259391e-01   
7.3690884894549e-01   
6.8852168077120e-01   
6.3685339445322e-01   
5.8215021256935e-01   
5.2467282046292e-01   
4.6469512391964e-01   
4.0250294385854e-01   
3.3839265425060e-01   
2.7266976975238e-01   
2.0564748978326e-01   
1.3764520598325e-01   
6.8986980163144e-02   
1.0816665110908e-16   ]},'coef',{[3.5826631552839e-03   
8.3231892962177e-03   
1.3031104991583e-02   
1.7677535257937e-02   
2.2239847550579e-02   
2.6696213967578e-02   
3.1025374934516e-02   
3.5206692201609e-02   
3.9220236729302e-02   
4.3046880709165e-02   
4.6668387718374e-02   
5.0067499237952e-02   
5.3228016731269e-02   
5.6134878759786e-02   
5.8774232718842e-02   
6.1133500831067e-02   
6.3201440073820e-02   
6.4968195750723e-02   
6.6425348449842e-02   
6.7565954163608e-02   
6.8384577378670e-02   
6.8877316977661e-02   
6.9041824829232e-02   ]});
data(10)= struct('node',{[9.9886640442007e-01   
9.9403196943209e-01   
9.8535408404801e-01   
9.7286438510669e-01   
9.5661095524281e-01   
9.3665661894488e-01   
9.1307855665579e-01   
8.8596797952361e-01   
8.5542976942995e-01   
8.2158207085934e-01   
7.8455583290040e-01   
7.4449430222607e-01   
7.0155246870682e-01   
6.5589646568544e-01   
6.0770292718495e-01   
5.5715830451465e-01   
5.0445814490746e-01   
4.4980633497404e-01   
3.9341431189757e-01   
3.3550024541944e-01   
2.7628819377953e-01   
2.1600723687604e-01   
1.5489058999815e-01   
9.3174701560086e-02   
3.1098338327189e-02   ]},'coef',{[2.9086225531552e-03   
6.7597991957456e-03   
1.0590548383651e-02   
1.4380822761485e-02   
1.8115560713489e-02   
2.1780243170125e-02   
2.5360673570013e-02   
2.8842993580535e-02   
3.2213728223579e-02   
3.5459835615146e-02   
3.8568756612588e-02   
4.1528463090147e-02   
4.4327504338803e-02   
4.6955051303949e-02   
4.9400938449466e-02   
5.1655703069581e-02   
5.3710621888997e-02   
5.5557744806213e-02   
5.7189925647728e-02   
5.8600849813222e-02   
5.9785058704266e-02   
6.0737970841770e-02   
6.1455899590316e-02   
6.1936067420683e-02   
6.2176616655347e-02   ]});
%% 積分變數的線性變換
[email protected](x)(((b+a)+(b-a)*x)/2);
%% 計算
area=0;
noden=data(n).node;
coefn=data(n).coef;
k=n;
n=n*h;
if n/2-floor(n/2)==1/2 %當階數為奇數   
   for i=1:floor(n/2)
      area=area+coefn(i)*(func(y(noden(i)))+func(y(-1*noden(i))));
   end
   i=floor(n/2)+1;
   area=area+coefn(i)*func(y(noden(i)));
elseif n/2-floor(n/2)==0 %當階數為偶數  
   for i=1:n/2
      area=area+coefn(i)*(func(y(noden(i)))+func(y(-1*noden(i))));
   end
end
area=area*(b-a)/2;
%% 誤差:前向（n-1階）差商作為事後誤差估計
if withError
    area_=gaussLegendreIntegral(func,a,b,k-1,false);
    gErrer=abs(area-area_);
end
end

輸出函式：

function str=printArray(a)
str=sprintf('%.13e\t',a);

4 Gauss積分的積分系數計算

4.1 測試Gauss積分的積分系數計算的函式

clc,clear
% GaussLegendre積分：(積分節點    積分系數)
% disp('GaussLegendre積分：(積分節點+-node    積分系數coef)')
% for n=1:16
%     sprintf('************階數 %d',n)
%     [x, w] = GaussLegendre_2(n);
%     k=floor(n/2);
%     if k~=n/2, spac=1:k+1;else, spac=1:n/2; end
%     node=-x(spac),coef=w(spac)
% end
%% GaussLaguerre積分：(積分節點    積分系數)
% disp('GaussLaguerre積分：(積分節點    積分系數)')
% alpha=1e-30
% for n=1:16
%     sprintf('************階數 %d',n)
%     [node,coef]=GaussLaguerre_2(n,alpha)
% end
%% GaussHermite 積分：(積分節點    積分系數)
% disp('GaussHermite積分：(積分節點    積分系數)')
% for n=1:20
%     sprintf('************階數 %d',n)
%     [x, w] =GaussHermite_2(n);
%     k=floor(n/2);
%     if k~=n/2, spac=1:k+1;else, spac=1:n/2; end
%     node=-x(spac),coef=w(spac)
% end
%%
printGaussIntData(1,5); 

printGaussIntData()函式（列印Gauss積分資料表）：

function printGaussIntData(type,h)
%% 列印Gauss積分資料表
str='';
for n=h:h:h*10
    if type==2
        [x, w] =GaussLaguerre_2(n,1e-100);
    else
        if type==1
            [x, w] =GaussLegendre_2(n);
        elseif type==3
            [x, w] =GaussHermite_2(n);
        end
        k=floor(n/2);
        if k~=n/2, spac=1:k+1;else, spac=1:n/2; end
        x=-x(spac); w=w(spac);
    end
    %data()= struct('node',{[]},'coef',{[]});sprintf('%d',n)
    str=[str,'data(',num2str((n/h)),...
        ')= struct(#&#node#&#,{[',printMatrix(x,'%.13e',false),...
        ']},#&#coef#&#,{[',printMatrix(w,'%.13e',false),']});',sprintf('\n')];
end,str
end

輸出矩陣的函式：

function s=printMatrix(m,style,with_endl)
if nargin<2,  style='%.13e'; end
if nargin<3,  with_endl=true; end
[dim1,dim2]=size(m);
s='';
for i=1:dim1
    s=[s,sprintf([style,'   '],m(i,:)),sprintf('\n')];
end
if ~with_endl
    s(end)='';
end
end

4.2 Gauss-Hermite積分系數求解函式

function [x, w] = GaussHermite_2(n)
%name：鄧能財     Date: 2013/12/23
i = 1:n-1; 
a = sqrt(i/2); 
CM = diag(a,1) + diag(a,-1);
[V L] = eig(CM); 
[x ind] = sort(diag(L)); 
V = V(:,ind)'; 
w = sqrt(pi) * V(:,1).^2;
end

4.3 Gauss-Laguerre積分系數求解函式

function [x, w] = GaussLaguerre_2(n, alpha)
%name：鄧能財     Date: 2013/12/23
i = 1:n; 
a = (2*i-1) + alpha; 
b = sqrt( i(1:n-1) .* ((1:n-1) + alpha) ); 
CM = diag(a) + diag(b,1) + diag(b,-1);
[V L] = eig(CM); 
[x ind] = sort(diag(L)); 
V = V(:,ind)'; 
w = gamma(alpha+1) .* V(:,1).^2;
end

4.4 Gauss-Legendre積分系數求解函式

function [x, w] = GaussLegendre_2(n)
% http://www.mathworks.cn/matlabcentral/fileexchange/8067-gauss-laguerre
%對應《數值分析方法 奚梅成的Gauss-Legendre積分》
%name：鄧能財     Date: 2013/12/23
i = 1:n-1; 
a = i./sqrt(4*i.^2-1); 
CM = diag(a,1) + diag(a,-1);
[V L] = eig(CM); 
[x ind] = sort(diag(L)); 
V = V(:,ind)'; 
w = 2 * V(:,1).^2;
end

5 NewtonCotes積分

5.1 測試NewtonCotes積分函式newtonCotesIntegral()

clc,clear
format long
disp('各個階數的Newton-Cote`s積分')
% [email protected](x)((x+1).^(-2))
f =  @(x)(cos(x))
a = -1.57079632679490
b = 1.57079632679490
% [email protected](x)(exp(-x.^2))
% [email protected](x)(exp(-x))
% [email protected](x)(exp(-x.^4))
% [email protected](x)((1/sqrt(2*pi))*exp(-x.^2/2)) %標準正太分佈
% [email protected](x)(sin(x)./x)
% [email protected](x)(x)
% a=0,b=1
area=[];nn=[];gErrer=[];
for n=2:8
    nn=[nn,n];
    [arean,gErrern]=newtonCotesIntegral(f,a,b,n,true);
    area=[area,arean];
    gErrer=[gErrer,gErrern];
end,nn,
area_str=printArray(area)
gErrer_str=printArray(gErrer)
%% 系統計算的積分精確值
area_system=quad(f,a,b,1e-13)
d_my_system=abs(area-area_system);
log_darea=round(log10(gErrer)*10)
log_d_my_system=round(log10(d_my_system)*10)
% %% 測試誤差項
% a=0,b=.01
% n=3
% [area,gErrer]=newtonCotesIntegral(f,a,b,n,true)
% area_system=quad(f,a,b,1e-13)

輸出陣列的函式：

function str=printArray(a)
str=sprintf('%.13e\t',a);

執行結果：

各個階數的Newton-Cote`s積分
f = 
    @(x)(cos(x))
a =
  -1.570796326794900
b =
   1.570796326794900
nn =
     2     3     4     5     6     7     8
area_str =
1.5707963267949e+00	2.2672492052928e+00	2.2888179796358e+00	2.2609034070862e+00	2.2389648201115e+00	2.2178021611519e+00	2.1987231887715e+00	
gErrer_str =
3.1341507153454e+00	1.7494105356846e+00	6.8206453262057e-02	1.1112991494925e-01	1.0995339715508e-01	1.3352735099341e-01	1.5154966451227e-01	
area_system =
     2

5.2 newtonCotesIntegral()

function [area,gErrer]=newtonCotesIntegral(func,a,b,n,withError)
%Newton-Cotes積分
%理論依據：數值分析方法 奚梅成
%func 被積函式；
%[a,b] 積分割槽間
%n Newton-Cotes積分的階數
%%
%name：鄧能財     Date: 2013/12/21
%% 預設引數
if nargin<5 withError=false; gErrer=0; end
%% 查錯
maxn=8;
% assert(a~=inf && b~=inf,'積分割槽間不能包含inf!');
% assert(ismember(n,[1:maxn]),['積分精確度數量級必須在1~',num2str(maxn),'之間!']);
%% 資料：1,2,...,maxn階Newton-Cotes積分的積分系數
    %   (積分系數coef)
data= struct('coef',cell(1,maxn));
%data()= struct('coef',{[]});
data(1)= struct('coef',{[1]});
data(2)= struct('coef',{[ 1/2          1/2  ]});
data(3)= struct('coef',{[ 1/6          4/6        1/6  ]});
data(4)= struct('coef',{[ 1/8          3/8        3/8        1/8 ]});
data(5)= struct('coef',{[ 7/90        16/45       2/15       16/45      7/90 ]});
data(6)= struct('coef',{[ 19/288      25/96      25/144      25/144     25/96      19/288 ]});
data(7)= struct('coef',{[ 41/840      9/35       9/280       34/105     9/280       9/35      41/840 ]});
data(8)= struct('coef',{[ 751/17280 3577/17280 1323/17280  2989/17280 2989/17280  1323/17280  3577/17280  751/17280 ]});
%% 計算
area=0;
coefn=data(n).coef;
xi=a; h=(b-a)/n;
for i=1:n
    xi=xi+h;
    area=area+coefn(i)*func(xi);
end
area=area*(b-a);
%% 誤差:前向（n-1階）差商作為事後誤差估計
if withError
    area_=newtonCotesIntegral(func,a,b,n-1,false);
    gErrer=abs(area-area_);
    gErrer=gErrer*(10^(.1*(n+1))); %根據試驗值調整誤差估計
end
% %% 誤差:兩分段（n-1階）積分和與原積分之差作為事後誤差估計
% if withError
%     area_=(newtonCotesIntegral(func,a,(a+b)/2,n,false)+...
%            newtonCotesIntegral(func,(a+b)/2,b,n,false));
%     gErrer=abs(area-area_);
% end
end

6 Romberg積分

6.1 測試Romberg積分函式rombergIntegral()

clc,clear
format long
f =  @(x)(cos(x))
a = -1.57079632679490
b = 1.57079632679490
'各種精度的Romberg積分'
% [email protected](x)(exp(x))
% a=1,b=2
area_=[];pre_=[];k_=[];
for i=2:5
    pre=10^(-2*i)
    [area,k]=rombergIntegral(f,a,b,pre);
    pre_=[pre_,pre];area_=[area_,area];k_=[k_,k];
end,
k_
area_str=printArray(area_)
pre_str=printArray(pre_)
area_system=quadl(f,a,b)

輸出陣列的函式：

function str=printArray(a)
str=sprintf('%.13e\t',a);

執行結果：

f = 
    @(x)(cos(x))
a =
  -1.570796326794900
b =
   1.570796326794900
ans =
各種精度的Romberg積分
pre =
     1.000000000000000e-04
pre =
     1.000000000000000e-06
pre =
     1.000000000000000e-08
pre =
     1.000000000000000e-10
k_ =
     4     5     6     6
area_str =
2.0000055499797e+00	1.9999999945873e+00	2.0000000000013e+00	2.0000000000013e+00	
pre_str =
1.0000000000000e-04	1.0000000000000e-06	1.0000000000000e-08	1.0000000000000e-10	
area_system =
   1.999999977471133

6.2 rombergIntegral()

function [area,k]=rombergIntegral(func,a,b,prec)
%Romberg積分
%理論依據：數值分析方法 奚梅成
%func 函式；
%[a,b] 積分割槽間
%pre 精度
%程式特點：密切結合理論，將記憶體和運算量達到極小化
    %只記錄最後一排需用的數字
%%
%name：鄧能財     Date: 2013/12/20
%這幾乎是我編寫的第一個【除錯只出一個錯誤】的程式
%% 計算
%取初始分段數為1
h=b-a;
t=[];%記錄最後一排需用的數字
%Romberg積分的階數：k
    %即所儲存的T的長度
k=1;
%計算梯形積分T0(h)
t=(func(a)+func(b))/2;
while true
    t=[0,t];
    k=k+1;
    %h1用於T0(h)與T0(h/2)之間的遞推計算
    h1=0;
    for xi=a+h/2:h:b-h/2
        h1=h1+func(xi);
    end, h1=h1*h;
    %計算梯形積分T0(h/2)
    t(1)=(t(2)+h1)/2;
    h=h/2;
    %計算所得積分最精確值
    bit=1;%記錄2^(2k)
    for i=2:k
        bit=bit*4;
        t(i)=t(i-1)+(t(i-1)-t(i))/(bit-1);
    end
    %判斷是否達到精度要求
    if abs(t(end)-t(end-1))<prec
        break;
    end
end
area=t(end);
end

7 復化積分

7.1 測試復化積分函式complexSimpson()

clc,clear
% [email protected](x)(exp(x))
% a=1,b=2
% m2=f(2)
% m4=m2
f =   @(x)(cos(x))
a = -1.57079632679490
b = 1.57079632679490
'各種精度的復化梯形積分積分'
% m2=1
m4=1
area_=[];pre_=[];n_=[];
for i=2:5
    pre=10^(-2*i);
    [area,n]=complexSimpson(f,a,b,pre,m4);
    pre_=[pre_,pre];area_=[area_,area];n_=[n_,n];
end,
n_
area_str=printArray(area_)
pre_str=printArray(pre_)
area_system=quadl(f,a,b)
% pre=1e-6
% ss=complexTrangleInt(f,a,b,pre,m2)
% ss=complexSimpson(f,a,b,pre,m4)

輸出陣列的函式：

function str=printArray(a)
str=sprintf('%.13e\t',a);

執行結果：

f = 
    @(x)(cos(x))
a =
  -1.570796326794900
b =
   1.570796326794900
ans =
各種精度的復化梯形積分積分
m4 =
     1
n_ =
     3     9    28    89
area_str =
2.0008631896735e+00	2.0000103477058e+00	2.0000001100950e+00	2.0000000010782e+00	
pre_str =
1.0000000000000e-04	1.0000000000000e-06	1.0000000000000e-08	1.0000000000000e-10	
area_system =
   1.999999977471133

7.2 complexSimpson()

function [ss,n]=complexSimpson(func,a,b,pre,m4)
%func函式
%a,b積分割槽間
%pre精度
%m4區間上4階導的最大值
%name：鄧能財     Date: 2013/12/23
%%%%%%%%%%%%%%%%%%%%%%：%【表需在尾部新增分號‘；’
n=floor((sqrt((b-a)^5)*m4/(2880*pre))^(1/4))+1%【
h=(b-a)/(2*n);
ss=func(a)+4*func(a+h)+func(b);
ss1=0; %n為奇數的項
ss2=0; %n為偶數的項
for xi=a+2*h:2*h:a+2*(n-1)*h
    ss1=ss1+func(xi+h);
    ss2=ss2+func(xi);
end
ss1,ss2
ss=ss+4*ss1+2*ss2;
ss=ss*(h/3);
end

7.3 complexTrangleInt()

function [ss,n]=complexTrangleInt(func,a,b,pre,m2)
%func函式
%a,b積分割槽間
%pre精度
%m2區間上二階導的最大值
%name：鄧能財     Date: 2013/12/23
%%%%%%%%%%%%%%%%%%%%%%：%【表需在尾部新增分號‘；’
n=floor((sqrt((b-a)^3)*m2/(12*pre))^.5)+1;%【
h=(b-a)/n;
ss=.5*(func(a)+func(b));
for xi=a+h:h:a+(n-1)*h
    ss=ss+func(xi);
end
ss=ss*h;
end

8 自適應運算積分

8.1 測試自適應運算積分函式selfAdaptIntegral()

clc,clear
format long
% [email protected](x)(exp(x))
% a=1,b=10
% pre=1e-13
f =  @(x)(cos(x))
a = -1.57079632679490
b = 1.57079632679490
% ss=complexTrangleInt(f,a,b,pre,m2)
% [area,deep]=selfAdaptIntegral(f,a,b,pre)
'各種精度的自適應Simpson積分'
area_=[];pre_=[];deep_=[];
for i=2:5
    pre=10^(-2*i);
    [area,deep]=selfAdaptIntegral(f,a,b,pre);
    pre_=[pre_,pre];area_=[area_,area];deep_=[deep_,deep];
end,
area_str=printArray(area_)
pre_str=printArray(pre_)
deep_
area_system=quadl(f,a,b)

輸出陣列的函式：

function str=printArray(a)
str=sprintf('%.13e\t',a);

執行結果：

f = 
    @(x)(cos(x))
a =
  -1.570796326794900
b =
   1.570796326794900
ans =
各種精度的自適應Simpson積分
area_str =
2.0000526243412e+00	2.0000002039922e+00	2.0000000017060e+00	2.0000000000498e+00	
pre_str =
1.0000000000000e-04	1.0000000000000e-06	1.0000000000000e-08	1.0000000000000e-10	
deep_ =
     1     3     5     6
area_system =
   1.999999977471133

8.2 selfAdaptIntegral()

function [area,deep]=selfAdaptIntegral(func,a,b,prec)
%Simpson公式的自適應積分
%理論依據：數值分析方法 奚梅成
%name：鄧能財     Date: 2013/12/23
%func 函式；
%[a,b] 積分割槽間
%pre 精度
%演算法：遞迴
%程式特點：密切結合理論，將記憶體和運算量達到極小化
%% 計算
%%disp('##Start')
%取分段次數為n
%n=3時，Simpson積分系數序列為: 1424241
			%，已經能夠發揮Simpson公式的作用，且分段數最少
n=3;
h=(b-a)/n;
%計算梯形積分T_n
s1=(func(a)+func(b))/2;
for xi=a+h:h:b-h
   s1=s1+func(xi);
end, s1=s1*h;
%H_n(f)
h1=0; 
for xi=a+h/2:h:b-h/2
   h1=h1+func(xi);
end, h1=h1*h;
%計算梯形積分T_2n
s2=(s1+h1)/2; 
%計算Simpson積分S_n
s1=(4*s2-s1)/3;
%H_2n(f)
h=h/2; h2=0; 
for xi=a+h/2:h:b-h/2
   h2=h2+func(xi);
end, h2=h2*h;
%計算Simpson積分S_2n
s2=s1/2+(4*h2-h1)/6;
%% 判斷是否達到要求精度
deep=1;
if(abs(s1-s2)<15*prec)
   area=s2;
   %disp('##End')
   return;
else %不滿足精度，進一步細分
    [area1,deep1]=selfAdaptIntegral(func,a,(a+b)/2,prec/2);
   	[area2,deep2]=selfAdaptIntegral(func,(a+b)/2,b,prec/2);
    area=area1+area2;
    deep=deep+max(deep1,deep2);
   %disp('##End')
   return;
end
end

9 積分系數——二

9.1 Gauss-Legendre積分系數

格式：序號【 積分節點 ******積分系數
1【
 0 ********* 2
2【
 ±0.5773502692  ********* 1
3【
 ±0.7745966692
0
 *********
0.5555555556
0.8888888889
4【
±0.8611363116
±0.3399810436
*********
0.3478548451
0.6521451549
5【
 ±0.9061798459
±0.5384693101
0
 *********
0.2369268851
0.4786286705
0.5688888889
6【
±0.9324695142
±0.6612093865
±0.2386191861
 *********
0.1713244924
0.3607615730
0.4679139346
7【
±0.9491079123
±0.7415311856
±0.4058451514
0
 *********
0.1294849662
0.2797053915
0.3818300505
0.4179591837
8【
±0.9602898565
±0.7966664774
±0.5255324099
±0.1834346425
 *********
0.1012285363
0.2223810345
0.3137066459
0.3626837834

9.2 newton cotes積分系數

表6—1
【】表錯誤
n  Ci(n-1)
 
2  1/2          1/2  
3  1/6          4/6        1/6  
4  1/8          3/8        3/8        1/8 
5  7/90        16/45       2/15       16/45      7/90 
6  19/288      25/96      25/144      25/144     25/96      19/288 
7  41/840      9/35       9/280       34/105     9/280       9/35      41/840 
8  751/17280 3577/17280 1323/17280  2989/17280 2989/17280  1323/17280  3577/17280  751/17280 
9  989/28350 5888/28350 -928/28350 10496/28350 -4540/28350 10496/28350 -928/28350 【5888/28350】 989/28350 

end