|
楼主 |
发表于 2012-1-4 09:40
|
显示全部楼层
没想到给大家添了这么多麻烦!
Gusssfun —— 进行高斯积分的函数
Guass_Point —— 区间[a,b]内进行高斯积分的积分点
grule —— 区间[-1,1]的高斯积分点
G_I —— 进行高斯积分点调用格式- function Y = Guassfun(u)
- Y = [1 + exp(-u)*sin(4*u);
- sin(sqrt(u))];
- % Y=[1 + u.^2 - 3 * u.^3 + 4 * u.^5, 1 + u.^2 - 3 * u.^3 ;
- % 1 + u.^2 - 3 * u.^3, u.^2 - 3 * u.^3 + 4 * u.^5;
- % ];
- end
复制代码- function [bp,wf] = grule(n)
- % -------------------------------------------------------------------------
- % This code is obtained by "Google" %
- % -------------------------------------------------------------------------
- % [bp,wf]=grule(n)
- % This function computes Gauss base points and weight factors
- % using the algorithm given by Davis and Rabinowitz in
- % 'Methods of Numerical Integration', page 365, Academic Press, 1975.
- % Input: n - the number of Guass integraltion points
- % Output: bp - the value of Guass integration points
- % wf - the value of Guass integration weights
- bp = zeros(n,1); wf = zeros(n,1);
- iter = 2; m = fix((n+1)/2); e1 = n*(n+1);
- mm = 4*m-1; t = (pi/(4*n+2))*(3:4:mm); nn = (1-(1-1/n)/(8*n*n));
- xo = nn*cos(t);
- for j = 1:iter
- pkm1 = 1; pk = xo;
- for k = 2:n
- t1 = xo.*pk;
- pkp1 = t1-pkm1-(t1-pkm1)/k+t1;
- pkm1 = pk; pk = pkp1;
- end
- den = 1.-xo.*xo; d1 = n*(pkm1-xo.*pk); dpn = d1./den;
- d2pn = (2.*xo.*dpn-e1.*pk)./den;
- d3pn = (4*xo.*d2pn+(2-e1).*dpn)./den;
- d4pn = (6*xo.*d3pn+(6-e1).*d2pn)./den;
- u = pk./dpn; v = d2pn./dpn;
- h = -u.*(1+(.5*u).*(v+u.*(v.*v-u.*d3pn./(3*dpn))));
- p = pk+h.*(dpn+(.5*h).*(d2pn+(h/3).*(d3pn+.25*h.*d4pn)));
- dp = dpn+h.*(d2pn+(.5*h).*(d3pn+h.*d4pn/3));
- h = h-p./dp; xo = xo+h;
- end
- bp = -xo-h;
- fx = d1-h.*e1.*(pk+(h/2).*(dpn+(h/3).*(d2pn+(h/4).*(d3pn+(.2*h).*d4pn))));
- wf = 2*(1-bp.^2)./(fx.*fx);
- if (m+m) > n, bp(m)=0; end
- if ~((m+m) == n), m=m-1; end
- jj = 1:m; n1j = (n+1-jj);
- bp(n1j) = -bp(jj); wf(n1j) = wf(jj);
- % -------------------------------------------------------------------------
- end
复制代码- function varginout = Guass_Point(n,a,b)
- % -------------------------------------------------------------------------
- %
- % Guass_Int.m
- %
- % Calculate the value of Guass integration points and weights for arbitrary interzone [a,b]
- %
- % Input : n - the number of Guass integraltion points
- % a - the low boundary of intergral interzone
- % b - the up boundary of intergral interzone
- % Output: varginout{1} = bp - the value of Guass integration points
- % varginout{2} = wf - the value of Guass integration weights
- %
- % Written by Yan - 21/3/2011
- % Contact: yanyongju@ase.buaa.edu.cn
- %
- % -------------------------------------------------------------------------
- % -------------------------------------------------------------------------
- % The Guass integral function can only used within the interzone [-1,1],
- % If we want to get the integration from integral boundary a to b, as [a,b]
- % we can tranform the [a,b] to [-1,1] with the equation
- % a + b b - a
- % s = ------- + ------- * t
- % 2 2
- % where t is in the interzone [-1,1], while s in [a,b]
- % so we can get
- % / a + b b - a \ b - a
- % f(s)ds = f| ------- + ------- * t | -------dt
- % \ 2 2 / 2
- % -------------------------------------------------------------------------
- if nargin == 1
- a = -1;
- b = 1;
- end
- [bp,wf] = grule(n);
- sp = size(bp);
- gpoint = (a+b)/2*ones(sp) + (b-a)/2 * bp;
- wpoint = (b-a)/2 * wf;
- varginout = {gpoint, wpoint};
- % -------------------------------------------------------------------------
- end
复制代码- function val = G_I(a,b,N)
- % a = 0; b = 1; N = 10;
- %
- % -------------------------------------------------------------------------
- %
- % Guass_Int.m
- %
- % Calculate the value of @Guassfun by Guass intergraltion
- %
- % Input : fun - the function of Guass integraltion
- % N - the number of Guass integraltion points
- % a - the low boundary of intergral interzone
- % b - the up boundary of intergral interzone
- % Output: val - the value of Guass integraltion
- %
- % Written by Yan - 21/3/2011
- % Contact: yanyongju@ase.buaa.edu.cn
- %
- % -------------------------------------------------------------------------
- par = Guass_Point(N,a,b);
- point = par{1}; weight = par{2};
- % -------------------------------------------------------------------------
- % func = Guassfun(point);
- % % func = feval('Guassfun',point);
- %
- % [m,n] = size(func);
- %
- % for i = 1:m
- % for j = 1:n
- % val(i,j) = dot(func{i,j},weight);
- % end
- % end
- % -------------------------------------------------------------------------
- val = 0.0;
- for k = 1:N
- func = Guassfun(point(k));
- val = val + func * weight(k);
- end
- % -------------------------------------------------------------------------
- end
复制代码 |
评分
-
1
查看全部评分
-
|