Function for initializing the matrices 'tt' and 'w' for Gaussian Integration


 % Function for initializing the matrices 'tt' and 'w' for Gaussian  
 % Integration; with n = 5, 10, 20, 24, 32, 40 or 48  
 % SAMPLE CODE FOR IMPLEMENTING DOUBLE INTEGRATION  
 % gaussian_points = 10 ; % Define the number of points to be used  
 %  
 % a = 1; % Lower integration limit for y  
 % b = 2; % Upper integration limit for y  
 %  
 % c = 3; % Lower integration limit for x  
 % d = 4; % Upper integration limit for x  
 %  
 % sum = 0;  
 %  
 % for ii = 1:gaussian_points  
 %  
 % yy = ((d - c)*tt(ii) + d + c)/2;  
 %  
 % inner_sum = 0;  
 %  
 % for jj = 1:gaussian_points  
 % xx = ((b - a)*tt(jj) + b + a)/2;  
 % Fm1 = xx + yy*yy; % I N T E G R A N D  
 % inner_sum = inner_sum + 0.5*(b - a)*Fm1*w(jj);  
 % end  
 %  
 % sum = sum + 0.5*(d - c)*inner_sum*w(ii);  
 %  
 % end  
 % -------------------------------------------------------------------------  
 function [tt w] = func_gaussian_integration_weights (gaussian_points)  
 switch (gaussian_points)  
 case 5  
 tt = [-0.906179845938664  
 -0.538469310105683  
 0.0  
 0.538469310105683  
 0.906179845938664]; % Points where function will be evaluated  
 w = [0.236926885056189  
 0.478628670499366  
 0.568888888888889  
 0.478628670499366  
 0.236926885056189]; % Constants to be multiplied  
 case 10  
 tt = [-0.973906528517172  
 -0.865063366688985  
 -0.679409568299024  
 -0.433395394129247  
 -0.148874338981631  
 0.148874338981631  
 0.433395394129247  
 0.679409568299024  
 0.865063366688985  
 0.973906528517172]; % Points where function will be evaluated  
 w = [0.066671344308688  
 0.149451349150581  
 0.219086362515982  
 0.269266719309996  
 0.295524224714753  
 0.295524224714753  
 0.269266719309996  
 0.219086362515982  
 0.149451349150581  
 0.066671344308688]; % Constants to be multiplied  
 case 20  
 tt = [-0.993128599185094924786  
 -0.963971927277913791268  
 -0.912234428251325905868  
 -0.839116971822218823395  
 -0.746331906460150792614  
 -0.636053680726515025453  
 -0.510867001950827098004  
 -0.373706088715419560673  
 -0.227785851141645078080  
 -0.076526521133497333755  
 0.076526521133497333755  
 0.227785851141645078080  
 0.373706088715419560673  
 0.510867001950827098004  
 0.636053680726515025453  
 0.746331906460150792614  
 0.839116971822218823395  
 0.912234428251325905868  
 0.963971927277913791268  
 0.993128599185094924786]; % Points where function will be evaluated  
 w = [0.017614007139152118312  
 0.040601429800386941331  
 0.062672048334109063570  
 0.083276741576704748725  
 0.101930119817240435037  
 0.118194531961518417312  
 0.131688638449176626898  
 0.142096109318382051329  
 0.149172986472603746788  
 0.152753387130725850698  
 0.152753387130725850698  
 0.149172986472603746788  
 0.142096109318382051329  
 0.131688638449176626898  
 0.118194531961518417312  
 0.101930119817240435037  
 0.083276741576704748725  
 0.062672048334109063570  
 0.040601429800386941331  
 0.017614007139152118312]; % Constants to be multiplied  
 case 24  
 tt = zeros (24, 1);  
 tt(13:24) = [0.064056892862605626085  
 0.191118867473616309159  
 0.315042679696163374387  
 0.433793507626045138487  
 0.545421471388839535658  
 0.648093651936975569252  
 0.740124191578554364244  
 0.820001985973902921954  
 0.866415527004401034213  
 0.938274552002732758524  
 0.974728555971309498198  
 0.995187219997021360180];  
 tt(1:12) = -1*flipud (tt (13:24)); % Points where function will be evaluated  
 w = zeros (24, 1);  
 w(13:24) = [0.127938195346752156974  
 0.125837456346828296121  
 0.121670472927803391204  
 0.115505668053725601353  
 0.107444270115965634783  
 0.097618652104113888270  
 0.086190161531953275917  
 0.073346481411080305734  
 0.059298584915436780746  
 0.044277438817419806169  
 0.028531388628933663181  
 0.012341229799987199547];  
 w(1:12) = flipud (w (13:24));  
 case 32  
 tt = zeros (32, 1);  
 tt(17:32) = [0.048307665687738316235  
 0.144471961582796493485  
 0.239287362252137074545  
 0.331868602282127649780  
 0.421351276130635345364  
 0.506899908932229390024  
 0.587715757240762329041  
 0.663044266930215200975  
 0.732182118740289680387  
 0.794483795967942406963  
 0.849367613732569970134  
 0.896321155766052123965  
 0.934906075937739689171  
 0.964762255587506430774  
 0.985611511545268335400  
 0.997263861849481563545];  
 tt(1:16) = -1*flipud (tt (17:32)); % Points where function will be evaluated  
 w = zeros (32, 1);  
 w(17:32) = [0.096540088514727800567  
 0.095638720079274859419  
 0.093844399080804565639  
 0.091173878695763884713  
 0.087652093004403811143  
 0.083311924226946755222  
 0.078193895787070306472  
 0.072345794108848506225  
 0.065822222776361846838  
 0.058684093478535547145  
 0.050998059262376176196  
 0.042835898022226680657  
 0.034273862913021433103  
 0.025392065309262059456  
 0.016274394730905670605  
 0.007018610009470096600];  
 w(1:16) = flipud (w (17:32));  
 case 40  
 tt = zeros (40, 1);  
 tt(21:40) = [0.038772417506050821933  
 0.116084070675255208483  
 0.192697580701371099716  
 0.268152185007253681141  
 0.341994090825758473007  
 0.413779204371605001525  
 0.483075801686178712909  
 0.549467125095128202076  
 0.612553889667980237953  
 0.671956684614179548379  
 0.727318255189927103281  
 0.778305651426519387965  
 0.824612230833311663196  
 0.865959503212259503281  
 0.902098806968874296728  
 0.932812808278676533361  
 0.957916819213791655805  
 0.977259949983774262663  
 0.990726238699457006453  
 0.998237709710559200350];  
 tt(1:20) = -1*flipud (tt (21:40)); % Points where function will be evaluated  
 w = zeros (40, 1);  
 w(21:40) = [0.077505947978424811264  
 0.077039818164247965588  
 0.076110361900626242372  
 0.074723169057968264200  
 0.072886582395804059061  
 0.070611647391286779695  
 0.067912045815233903826  
 0.064804013456601038075  
 0.061306242492928939167  
 0.057439769099391551367  
 0.053227846983936824355  
 0.048695807635072232061  
 0.043870908185673271992  
 0.038782167974472017640  
 0.033460195282547847393  
 0.027937006980023401098  
 0.022245849194166957262  
 0.016421058381907888713  
 0.010498284531152813615  
 0.004521277098533191258];  
 w(1:20) = flipud (w (21:40));  
 case 48  
 tt = zeros (48, 1);  
 tt(25:48) = [0.032380170962869362033  
 0.097004699209462698930  
 0.161222356068891718056  
 0.224763790394689061225  
 0.287362487355455576736  
 0.348755886292160738160  
 0.408686481990716729916  
 0.466902904750958404545  
 0.523160974722233033678  
 0.577224726083972703818  
 0.628867396776513623995  
 0.677872379632663905212  
 0.724034130923814654674  
 0.767159032515740339254  
 0.807066204029442627083  
 0.843588261624393530711  
 0.876572020274247885906  
 0.905879136715569672822  
 0.931386690706554333114  
 0.952987703160430860723  
 0.970591592546247250461  
 0.984124583722826857745  
 0.993530172266350757548  
 0.998771007252426118601];  
 tt(1:24) = -1*flipud (tt (25:48)); % Points where function will be evaluated  
 w = zeros (48, 1);  
 w(25:48) = [0.064737696812683922503  
 0.064466164435950082207  
 0.063924238584648186624  
 0.063114192286254025657  
 0.062039423159892663904  
 0.06070443916589388053  
 0.059114839698395635746  
 0.057277292100403215705  
 0.055199503699984162868  
 0.052890189485193667096  
 0.050359035553854474958  
 0.047616658492490474826  
 0.044674560856694280419  
 0.041545082943464749214  
 0.038241351065830706317  
 0.034777222564770438893  
 0.031167227832798088902  
 0.027426509708356948200  
 0.023570760839324379141  
 0.019616160457355527814  
 0.015579315722943848728  
 0.011477234579234539490  
 0.007327553901276262102  
 0.003153346052305838633];  
 w(1:24) = flipud (w (25:48));  
 otherwise  
 error('Parameter "gaussian_points" can only take values 5, 10, 20, 24, 32, 40 or 48');  
 end  
 end