From 882a25b3360e18554912947e3d2f8a1c7a11fb61 Mon Sep 17 00:00:00 2001 From: treecity Date: Wed, 13 Apr 2011 18:14:08 +0000 Subject: [PATCH] =?utf8?q?[gGF]=20Parameter=20aufgesplittet=20[src]=20test?= =?utf8?q?=5Fsolve=20hinzugef=C3=BCgt.=20funktioniert=20aber=20noch=20nich?= =?utf8?q?t?= MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@8 26120e32-c555-405d-b3e1-1f783fb42516 --- src/F_par.m | 17 +++++++++++------ src/G00.m | 23 ++++++++++++----------- src/g0.m | 3 ++- src/test_solve.m | 37 +++++++++++++++++++++++++++++++++++++ 4 files changed, 62 insertions(+), 18 deletions(-) create mode 100644 src/test_solve.m diff --git a/src/F_par.m b/src/F_par.m index cb7d43b..a1794aa 100644 --- a/src/F_par.m +++ b/src/F_par.m @@ -1,11 +1,16 @@ -function sol = F_par(x, y, delta) +function sol = F_par(x1,x2,y1,y2,d1,d2,d3) - xyd=x-y-delta(1:2); +% xyd=x-y-delta(1:2); + + sol = (x1-y1-d1)*(x2-y2-d2)*G00(-0.5,x1,x2,y1+d1,y2+d2,d3)... + -(x1-y1-d1)*g0(0.5,x1,x1+d1,sqrt((x2-y2-d2)^2+d3^2))... + -(x2-y2-d2)*g0(0.5,x2,x2+d2,sqrt((x1-y1-d1)^2+d3^2))... + +1/3*((x1-y1-d1)^2+(x2-y2-d2)^2+d3^2)^(1.5); - sol = xyd(1)*xyd(2)*G00(-0.5,x,y+delta(1:2),delta(3)... - -xyd(1)*g0(0.5,x(1),y(1)+delta(1),norm([xyd(2) delta(3)]))... - -xyd(2)*g0(0.5,x(2),y(2)+delta(2),norm([xyd(1) delta(3)]))... - +1/3*(xyd(1)^2+xyd(2)^2+delta(3)^2)^(1.5)); +% sol = xyd(1)*xyd(2)*G00(-0.5,x,y+delta(1:2),delta(3))... +% -xyd(1)*g0(0.5,x(1),y(1)+delta(1),norm([xyd(2) delta(3)]))... +% -xyd(2)*g0(0.5,x(2),y(2)+delta(2),norm([xyd(1) delta(3)]))... +% +1/3*(xyd(1)^2+xyd(2)^2+delta(3)^2)^(1.5); end diff --git a/src/G00.m b/src/G00.m index dd0f40d..f181612 100644 --- a/src/G00.m +++ b/src/G00.m @@ -1,18 +1,19 @@ -function sol = G00(p,y,x,lambda) +function sol = G00(p,y1,y2,x1,x2,l) - sol = nan; - yx = y-x; +% sol = nan; if(p==-1.5) - if(lambda==0) - sol = -norm(yx)/prod(yx); + if(l==0) + sol = -sqrt((y1-x1)^2+(y2-x2)^2)/((y1-x1)*(y2-x2)); else - sol = sign(prod(yx))/(2*abs(lambda))*... - acos(-2*prod(yx)^2/prod((yx).^2+lambda^2)+1); + sol = sign((y1-x1)*(y2-x2))/(2*abs(l)) ... + * acos(-2*(y1-x1)^2*(y2-x2)^2 ... + /(((y1-x1)^2+l^2)*((y2-x2)^2+l^2))+1); end elseif(p==-0.5) - sol = (2*p*lambda^2*G00(p-1,y,x,lambda)... - + sum( yx(1).*g0(p,y(2),x(2),norm([yx(1) lambda])))... - + sum( yx(2).*g0(p,y(1),x(1),norm([yx(2) lambda]))))/(2*p+2); + sol = 2*p*l^2*(G00(p-1,y1,y2,x1,x2,l) ... + + (y1-x1)*g0(p,y2,x2,sqrt((y1-x1)^2+l^2)) ... + + (y2-x2)*g0(p,y1,x1,sqrt((y2-x2)^2+l^2)))/(2*p+2); + else + error 'no case for p defined' end - end diff --git a/src/g0.m b/src/g0.m index c8e1632..d576073 100644 --- a/src/g0.m +++ b/src/g0.m @@ -1,6 +1,6 @@ function sol = g0(p, y, x, lambda) - sol = nan; + %sol = nan; yx = y-x; if(lambda~=0) if(p==0.5) @@ -23,6 +23,7 @@ function sol = g0(p, y, x, lambda) sol = (yx).*abs(yx).^(2*p)/(2*p+1); end end + end function sol = arsinh(x) %inline diff --git a/src/test_solve.m b/src/test_solve.m new file mode 100644 index 0000000..4aa6aed --- /dev/null +++ b/src/test_solve.m @@ -0,0 +1,37 @@ + +N = size(elements); + +A = zeros(N); + + + +% intF = @(f,k1,k2,l1,l2) ... +% f(k1,k2,l1,l2)-f(k1,k2,l1,0)-f(k1,k2,0,l2)+f(k1,k2,0,0)... +% -f(k1,0,l1,l2)+f(k1,0,l1,0)+f(k1,0,0,l2)-f(k1,0,0,0)... +% -f(0,k2,l1,l2)+f(0,k2,l1,0)+f(0,k2,0,l2)-f(0,l2,0,0)... +% +f(0,0,l1,l2)-f(0,0,l1,0)-f(0,0,0,l2)+f(0,0,0,0); +% +% intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3)),x(1),x(2),y(1),y(2)); +for j = 1:N + + + for k = j+1:N + + + ej = coordinates(elements(j,:)',:); + ek = coordinates(elements(k,:)',:); + d = ek(1,:) - ej(1,:) + + x = [sum(ej(2,:)-ej(1,:)) sum(ej(3,:)-ej(1,:))] + y = [sum(ek(2,:)-ek(1,:)) sum(ek(3,:)-ek(1,:))] + + A(j,k) = 1/ (4*pi) *... intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3)),x(1),x(2),y(1),y(2)); + F_par(x(1),x(2),y(1),y(2),d(1),d(2),d(3)); + end +end + +A +b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2)) + + +x = A\b; \ No newline at end of file -- 2.47.3