\r
if (rx = ry) {\r
if (rx = 2) {\r
- //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x1[0],x1[1],x2[0],x3[1],y1[0],y1[1],y2[0],y3[1],d[0],d[1],d[2]);\r
+ printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x1[0], x2[0], x1[1], x3[1], y1[0],\r
+ y2[0], y1[1], y3[1], d[0], d[1], d[2]);\r
//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
tmp = quadInt(F_par, x1[0], x2[0], x1[1], x3[1], y1[0],\r
y2[0], y1[1], y3[1], d[0], d[1], d[2]);\r
ek = coordinates(elements(k,:)',:);
d = ek(1,:) - ej(1,:);
% d = ones(1,3);
- A1(j,k) = intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
- ,ej(1,1),ej(2,1),ej(1,2), ej(3,2),ek(1,1), ek(2,1),ek(1,2), ek(3,2));
+ A1(j,k) = surfDoubleQuad(@(x1,x2,y1,y2)1/sqrt((x1-y1-d(1)).^2+(x2-y2-d(2)).^2+d(3).^2)...
+ ,ej(1,1),ej(2,1),ej(1,2), ej(3,2),ek(1,1), ek(2,1),ek(1,2), ek(3,2),4);
+
+% A1(j,k) = intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
+% ,ej(1,1),ej(2,1),ej(1,2), ej(3,2),ek(1,1), ek(2,1),ek(1,2), ek(3,2));
end
end
--- /dev/null
+function [nodes,weights] = gauss(n,varargin)
+
+% Compute Gauss-Legendre quadrature on compact interval [a,b]
+%
+% Usage: [NODES,WEIGHTS] = gauss(N [,A,B])
+%
+% Here, N is the length of the Gaussian quadrature rule, and the
+% optional scalars A and B determine the compact interval. By default,
+% there holds A = -1 and b = +1.
+%
+% The function returns an (N x 1) column vector WEIGHTS containing the
+% quadrature weights and an (1 x N) row vector NODES containing the
+% corresponding quadrature nodes.
+%
+% Example: Suppose we aim to approximate the integral
+%
+% int_a^b f dx
+%
+% Assume that F is a Matlab function, which takes a column vector of
+% evaluation points X and returns a column vector Y of function values,
+% where Y(j) = F(X(j)). Then, the numerical integration reads
+%
+% [NODES,WEIGHTS] = gauss(n,A,B);
+% INT = WEIGHTS*F(NODES);
+%
+% Author: Dirk Praetorius - 11.03.2010
+
+
+beta = (1:n-1)./sqrt((2*(1:n-1)).^2-1);
+A = diag(beta,-1)+diag(beta,1);
+[eigenvector,nodes] = eig(A);
+[nodes,idx] = sort(diag(nodes));
+weights = 2*eigenvector(1,idx).^2;
+
+if nargin >= 3
+ a = varargin{1};
+ b = varargin{2};
+ weights = 0.5*abs(b-a)*weights;
+ nodes = 0.5*( a+b + nodes*(b-a) );
+end
--- /dev/null
+
+mex mex_g0.cpp SLPrecangle.cpp
+mex mex_G00.cpp SLPrecangle.cpp
+mex mex_Fpar.cpp SLPrecangle.cpp
+
+mex build_A.cpp SLPrecangle.cpp
+
+% function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
+% %Vierfachintegral ueber f(x1,x2,y1,y2) auf [a b]x[c d]x[r t]x[u v] mit
+% %der Simpsonregel und w Auswertungspunkten
+% s = (b-a)/w;
+% in = 0;
+% for i = a:s:b-s
+% in = in + fx2(f,i,c,d,r,t,u,v,w) + ...
+% 4*fx2(f,(2*i+s)/2,c,d,r,t,u,v,w) + ...
+% fx2(f,i+s,c,d,r,t,u,v,w);
+% end
+% in = in *s/6;
+% end
+%
+% function in = fx2(f,x1,c,d,r,t,u,v,w)
+% s = (d-c)/w;
+% in = 0;
+% for i = c:s:d-s
+% in = in + fy1(f,x1,i,r,t,u,v,w) + ...
+% 4*fy1(f,x1,(2*i+s)/2,r,t,u,v,w) + ...
+% fy1(f,x1,i+s,r,t,u,v,w);
+% end
+% in = in *s/6;
+% end
+%
+% function in = fy1(f,x1,x2,r,t,u,v,w)
+% s = (t-r)/w;
+% in = 0;
+% for i = r:s:t-s
+% in = in + fy2(f,x1,x2,i,u,v,w) + ...
+% 4*fy2(f,x1,x2,(2*i+s)/2,u,v,w) + ...
+% fy2(f,x1,x2,i+s,u,v,w);
+% end
+% in = in *s/6;
+% end
+%
+% function in = fy2(f,x1,x2,y1,u,v,w)
+% s = (v-u)/w;
+% in = 0;
+% for i = u:s:v-s
+% in = in + f(x1,x2,y1,i) + ...
+% 4*f(x1,x2,y1,(2*i+s)/2) + ...
+% f(x1,x2,y1,i+s);
+% end
+% in = in *s/6;
+% end
+
+
function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
-%Vierfachintegral ueber f(x1,x2,y1,y2) auf [a b]x[c d]x[r t]x[u v] mit
-%der Simpsonregel und w Auswertungspunkten
- s = (b-a)/w;
- in = 0;
- for i = a:s:b-s
- in = in + fx2(f,i,c,d,r,t,u,v,w) + ...
- 4*fx2(f,(2*i+s)/2,c,d,r,t,u,v,w) + ...
- fx2(f,i+s,c,d,r,t,u,v,w);
- end
- in = in *s/6;
-end
-function in = fx2(f,x1,c,d,r,t,u,v,w)
- s = (d-c)/w;
- in = 0;
- for i = c:s:d-s
- in = in + fy1(f,x1,i,r,t,u,v,w) + ...
- 4*fy1(f,x1,(2*i+s)/2,r,t,u,v,w) + ...
- fy1(f,x1,i+s,r,t,u,v,w);
- end
- in = in *s/6;
-end
+ [x1n x1w] = gauss(w,a,b);
+ [x2n x2w] = gauss(w,c,d);
+ [y1n y1w] = gauss(w,r,t);
+ [y2n y2w] = gauss(w,u,v);
-function in = fy1(f,x1,x2,r,t,u,v,w)
- s = (t-r)/w;
in = 0;
- for i = r:s:t-s
- in = in + fy2(f,x1,x2,i,u,v,w) + ...
- 4*fy2(f,x1,x2,(2*i+s)/2,u,v,w) + ...
- fy2(f,x1,x2,i+s,u,v,w);
+ for i=1:length(y2n)
+ for j=1:length(y1n)
+ for k=1:length(x2n)
+ for l = 1:length(x1n)
+
+ in = in + y2w(i) * y1w(j) * x2w(k) * x1w(l) * f(x1n(l),x2n(k),y1n(j),y2n(i));
+
+ end
+ end
+ end
end
- in = in *s/6;
-end
-function in = fy2(f,x1,x2,y1,u,v,w)
- s = (v-u)/w;
- in = 0;
- for i = u:s:v-s
- in = in + f(x1,x2,y1,i) + ...
- 4*f(x1,x2,y1,(2*i+s)/2) + ...
- f(x1,x2,y1,i+s);
- end
- in = in *s/6;
end
\ No newline at end of file
+% function in = surfQuad(f,a,b,c,d,v)\r
+% %Doppelintegral ueber f(x,y) auf [a b]x[c d] mit der Simpsonregel und v\r
+% %Auswertungspunkten\r
+% s = (b-a)/v;\r
+% in = 0;\r
+% for i = a:s:b-s\r
+% in = in + fy(f,i,c,d,v) + 4*fy(f,(2*i+s)/2,c,d,v) + fy(f,i+s,c,d,v);\r
+% end\r
+% in = in *s/6;\r
+% end\r
+% \r
+% function in = fy(f,x,c,d,v)\r
+% s = (d-c)/v;\r
+% in = 0;\r
+% for i = c:s:d-s\r
+% in = in + f(x,i) + 4*f(x,(2*i+s)/2) + f(x,i+s);\r
+% end\r
+% in = in *s/6;\r
+% end\r
+\r
+\r
function in = surfQuad(f,a,b,c,d,v)\r
-%Doppelintegral ueber f(x,y) auf [a b]x[c d] mit der Simpsonregel und v\r
-%Auswertungspunkten\r
- s = (b-a)/v;\r
+ [xn xw] = gauss(v,a,b);\r
+ [yn yw] = gauss(v,c,d);\r
+ \r
in = 0;\r
- for i = a:s:b-s\r
- in = in + fy(f,i,c,d,v) + 4*fy(f,(2*i+s)/2,c,d,v) + fy(f,i+s,c,d,v);\r
+ for i = 1:length(yn)\r
+ in = in + yw(i) * xw*f(xn,yn(i));\r
end\r
- in = in *s/6;\r
-end\r
\r
-function in = fy(f,x,c,d,v)\r
- s = (d-c)/v;\r
- in = 0;\r
- for i = c:s:d-s\r
- in = in + f(x,i) + 4*f(x,(2*i+s)/2) + f(x,i+s);\r
- end\r
- in = in *s/6;\r
end
\ No newline at end of file
done = 0;
%% g0 testen
p_g0 = [0.5 0 -0.5 -1 -1.5];
-
+for l= 1:5
for p = p_g0
- s = test_g0(p,1.4,2,1,4);
+ s = test_g0(p,1.4,l,1,4);
if(abs(s(1)-s(2))>eps)
+ disp g0
+ l
p
s
return;
end
end
-
+end
% x muss ausserhalb des Intervalls liegen
for p = p_g0
s = test_g0(p,1.4,0,2,4);
if(abs(s(1)-s(2))>eps)
+ disp g0
p
s
return;
for p = p_G00
s = test_G00(p,1.4,0.7,1,2,4,1,3);
if(abs(s(1)-s(2))>eps)
+ disp G00
p
s
return;
for p = p_G00
s = test_G00(p,1.1,0.3,0,2,4.5,1,2.5);
if(abs(s(1)-s(2))>eps)
+ disp G00
p
s
return;
s = test_Fpar(1,3,2,3,4,5,1,3,1,3,1);
if(abs(s(1)-s(2))>eps)
+ disp Fpar
p
s
return;
function [ret] = test_g0(p, x,l,a,b)
-quad_sol = quad(@(y)((x-y).^2+l.^2).^p,a,b);
-
-X = 0:0.2:5;
-Y = ((x-X).^2+l.^2).^p;
-
-plot(X,Y)
+% quad_sol = quad(@(y)((x-y).^2+l.^2).^p,a,b);
+[no we] = gauss(11,a,b);
+f = @(y)((x-y).^2+l.^2).^p;
+quad_sol = we*f(no);
+
+% X = 0:0.2:5;
+% Y = ((x-X).^2+l.^2).^p;
+%
+% plot(X,Y)
[my(1) mk(1)] = mex_g0(p,a,x,l);
[my(2) mk(2)] = mex_g0(p,b,x,l);
load exmpl_2DLShape\r
\r
\r
-[coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
-[coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
+% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
+% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
\r
+% build_A(coordinates,elements);\r
\r
-A = build_A(coordinates,elements);\r
+A = build_A2(coordinates,elements)\r
sum(sum (A-A'))\r
b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));\r
x = A\b\r
-xe = x'*A*x;\r
+xe = x'*A*x\r
\r
[coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
-A_fine = build_A(coordinates_fine,elements_fine);\r
+A_fine = build_A2(coordinates_fine,elements_fine)\r
sum(sum(A_fine-A_fine'))\r
b = sqrt(sum(quad_norm(coordinates_fine,elements_fine,'w').^2,2));\r
x_fine = A_fine\b;\r
projects / bacc.git / commitdiff