|
0
|
1 function int = quadg(fun,xlow,xhigh,tol,trace,p1,p2,p3,p4,p5,p6,p7,p8,p9) |
|
|
2 % |
|
|
3 %usage: int = quadg('Fun',xlow,xhigh) |
|
|
4 %or |
|
|
5 % int = quadg('Fun',xlow,xhigh,tol) |
|
|
6 %or |
|
|
7 % int = quadg('Fun',xlow,xhigh,tol,trace,p1,p2,....) |
|
|
8 % |
|
|
9 %This function works just like QUAD or QUAD8 but uses a Gaussian quadrature |
|
|
10 %integration scheme. Use this routine instead of QUAD or QUAD8: |
|
|
11 % if higher accuracy is desired (this works best if the function, |
|
|
12 % 'Fun', can be approximated by a power series) |
|
|
13 % or if many similar integrations are going to be done (I think less |
|
|
14 % function evaluations will typically be done, but the |
|
|
15 % integration points and the weights must be calculated. |
|
|
16 % These are saved between integrations so when QUADG |
|
|
17 % is called again, the points and weights are all ready |
|
|
18 % known.) |
|
|
19 % or if the function evaluations are time consuming. |
|
|
20 %Note that if there are discontinuities the integral should be broken up into separate |
|
|
21 %pieces. And if there are singularities, a more appropriate integration quadrature |
|
|
22 %should be used (such as the Gauss-Chebyshev). |
|
|
23 |
|
|
24 global b2 |
|
|
25 global w2 |
|
|
26 |
|
|
27 if ( exist('tol') != 1 ) |
|
|
28 tol=1e-3; |
|
|
29 elseif ( tol==[] ) |
|
|
30 tol=1e-3; |
|
|
31 endif |
|
|
32 if ( exist('trace') != 1 ) |
|
|
33 trace=0; |
|
|
34 elseif ( trace==[] ) |
|
|
35 trace=0; |
|
|
36 else |
|
|
37 trace=1; |
|
|
38 endif |
|
|
39 |
|
|
40 %setup string to call the function |
|
|
41 exec_string=['y=',fun,'(x']; |
|
|
42 num_parameters=nargin-5; |
|
|
43 for i=1:num_parameters |
|
|
44 exec_string=[exec_string,',p',int2str(i)]; |
|
|
45 endfor |
|
|
46 exec_string=[exec_string,');']; |
|
|
47 |
|
|
48 %setup mapping parameters |
|
|
49 jacob=(xhigh-xlow)/2; |
|
|
50 |
|
|
51 %generate the first two sets of integration points and weights |
|
|
52 if ( exist('b2') != 1 ) |
|
|
53 [b2,w2]=grule(2); |
|
|
54 endif |
|
|
55 |
|
|
56 x=(b2+1)*jacob+xlow; |
|
|
57 eval(exec_string); |
|
|
58 int_old=sum(w2.*y)*jacob; |
|
|
59 if ( trace==1 ) |
|
|
60 x_trace=x(:); |
|
|
61 y_trace=y(:); |
|
|
62 endif |
|
|
63 |
|
|
64 converge=0; |
|
|
65 for i=1:7 |
|
|
66 gnum=int2str(2^(i+1)); |
|
|
67 vname = ['b',gnum]; |
|
|
68 if ( exist(vname) == 0 ) |
|
|
69 estr = ['[b',gnum,',w',gnum,']=grule(',gnum,');']; |
|
|
70 eval(estr); |
|
|
71 estr = ['global b',gnum,' w',gnum,';']; |
|
|
72 eval(estr); |
|
|
73 endif |
|
|
74 estr = ['x=(b',gnum,'+1)*jacob+xlow;']; |
|
|
75 eval(estr); |
|
|
76 eval(exec_string); |
|
|
77 estr = ['int=sum(w',gnum,'.*y)*jacob;']; |
|
|
78 eval(estr); |
|
|
79 |
|
|
80 if ( trace==1 ) |
|
|
81 x_trace=[x_trace;x(:)]; |
|
|
82 y_trace=[y_trace;y(:)]; |
|
|
83 endif |
|
|
84 |
|
|
85 if ( abs(int_old-int) < abs(tol*int) ) |
|
|
86 converge=1; |
|
|
87 break; |
|
|
88 endif |
|
|
89 int_old=int; |
|
|
90 endfor |
|
|
91 |
|
|
92 if ( converge==0 ) |
|
|
93 disp('Integral did not converge--singularity likely') |
|
|
94 endif |
|
|
95 |
|
|
96 if ( trace==1 ) |
|
|
97 plot(x_trace,y_trace,'+') |
|
|
98 endif |
|
|
99 |
|
|
100 %gnum,i,length(x_trace) |
|
|
101 |
|
|
102 endfunction |
