From e2eb2a850b66e12dae5fee579ec9ec329e8d8e11 Mon Sep 17 00:00:00 2001 From: treecity Date: Mon, 16 May 2011 07:47:48 +0000 Subject: [PATCH] [src] Gauss-Quadratur zum Testen [src] test_functions zum Testen von g G und F [src] mexit zum Bauen der MEX Funktionen git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@25 26120e32-c555-405d-b3e1-1f783fb42516 --- src/build_A.cpp | 3 +- src/build_A2.m | 7 ++- src/gauss.m | 40 +++++++++++++++++ src/mexit.m | 7 +++ src/surfDoubleQuad.m | 99 ++++++++++++++++++++++++++----------------- src/surfQuad.m | 40 +++++++++++------ src/test_functions.m | 12 ++++-- src/test_g0.m | 15 ++++--- src/test_solveError.m | 11 ++--- 9 files changed, 165 insertions(+), 69 deletions(-) create mode 100644 src/gauss.m create mode 100644 src/mexit.m diff --git a/src/build_A.cpp b/src/build_A.cpp index 03c5054..9281525 100644 --- a/src/build_A.cpp +++ b/src/build_A.cpp @@ -114,7 +114,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (rx = ry) { if (rx = 2) { - //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]); + 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], + y2[0], y1[1], y3[1], d[0], d[1], d[2]); //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]); tmp = quadInt(F_par, x1[0], x2[0], x1[1], x3[1], y1[0], y2[0], y1[1], y3[1], d[0], d[1], d[2]); diff --git a/src/build_A2.m b/src/build_A2.m index 9caa83a..9ca4079 100644 --- a/src/build_A2.m +++ b/src/build_A2.m @@ -16,8 +16,11 @@ function A = build_A2(coordinates,elements) 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 diff --git a/src/gauss.m b/src/gauss.m new file mode 100644 index 0000000..795c59f --- /dev/null +++ b/src/gauss.m @@ -0,0 +1,40 @@ +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 diff --git a/src/mexit.m b/src/mexit.m new file mode 100644 index 0000000..4057572 --- /dev/null +++ b/src/mexit.m @@ -0,0 +1,7 @@ + +mex mex_g0.cpp SLPrecangle.cpp +mex mex_G00.cpp SLPrecangle.cpp +mex mex_Fpar.cpp SLPrecangle.cpp + +mex build_A.cpp SLPrecangle.cpp + diff --git a/src/surfDoubleQuad.m b/src/surfDoubleQuad.m index 49a9e3e..17fed8f 100644 --- a/src/surfDoubleQuad.m +++ b/src/surfDoubleQuad.m @@ -1,45 +1,68 @@ +% 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 diff --git a/src/surfQuad.m b/src/surfQuad.m index b089675..3111d7d 100644 --- a/src/surfQuad.m +++ b/src/surfQuad.m @@ -1,19 +1,31 @@ +% function in = surfQuad(f,a,b,c,d,v) +% %Doppelintegral ueber f(x,y) auf [a b]x[c d] mit der Simpsonregel und v +% %Auswertungspunkten +% s = (b-a)/v; +% in = 0; +% for i = a:s:b-s +% 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); +% end +% in = in *s/6; +% end +% +% function in = fy(f,x,c,d,v) +% s = (d-c)/v; +% in = 0; +% for i = c:s:d-s +% in = in + f(x,i) + 4*f(x,(2*i+s)/2) + f(x,i+s); +% end +% in = in *s/6; +% end + + function in = surfQuad(f,a,b,c,d,v) -%Doppelintegral ueber f(x,y) auf [a b]x[c d] mit der Simpsonregel und v -%Auswertungspunkten - s = (b-a)/v; + [xn xw] = gauss(v,a,b); + [yn yw] = gauss(v,c,d); + in = 0; - for i = a:s:b-s - 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); + for i = 1:length(yn) + in = in + yw(i) * xw*f(xn,yn(i)); end - in = in *s/6; -end -function in = fy(f,x,c,d,v) - s = (d-c)/v; - in = 0; - for i = c:s:d-s - in = in + f(x,i) + 4*f(x,(2*i+s)/2) + f(x,i+s); - end - in = in *s/6; end \ No newline at end of file diff --git a/src/test_functions.m b/src/test_functions.m index 00ba036..ce33b3a 100644 --- a/src/test_functions.m +++ b/src/test_functions.m @@ -3,20 +3,23 @@ function done=test_functions(eps) 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; @@ -29,6 +32,7 @@ p_G00 = [-0.5 -1.5]; 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; @@ -38,6 +42,7 @@ p_G00 = [-0.5 -1.5]; 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; @@ -49,6 +54,7 @@ p_G00 = [-0.5 -1.5]; 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; diff --git a/src/test_g0.m b/src/test_g0.m index 1bd8d15..89e2a02 100644 --- a/src/test_g0.m +++ b/src/test_g0.m @@ -1,12 +1,15 @@ 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); diff --git a/src/test_solveError.m b/src/test_solveError.m index 9a1c633..0049235 100644 --- a/src/test_solveError.m +++ b/src/test_solveError.m @@ -3,18 +3,19 @@ load exmpl_2DLShape -[coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1))); -[coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1))); +% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1))); +% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1))); +% build_A(coordinates,elements); -A = build_A(coordinates,elements); +A = build_A2(coordinates,elements) sum(sum (A-A')) b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2)); x = A\b -xe = x'*A*x; +xe = x'*A*x [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1))); -A_fine = build_A(coordinates_fine,elements_fine); +A_fine = build_A2(coordinates_fine,elements_fine) sum(sum(A_fine-A_fine')) b = sqrt(sum(quad_norm(coordinates_fine,elements_fine,'w').^2,2)); x_fine = A_fine\b; -- 2.47.3