From 2bcc59674c881d890682f9db26cd2168f14b6ad1 Mon Sep 17 00:00:00 2001 From: maria Date: Sat, 16 Apr 2011 09:29:42 +0000 Subject: [PATCH] [src] Funktionen laufen jetzt git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@10 26120e32-c555-405d-b3e1-1f783fb42516 --- src/F_par.m | 25 +++++++++++++++++++++---- src/G00.m | 13 +++++++++---- src/test_solve.m | 15 +++------------ 3 files changed, 33 insertions(+), 20 deletions(-) diff --git a/src/F_par.m b/src/F_par.m index 95b7842..9a57744 100644 --- a/src/F_par.m +++ b/src/F_par.m @@ -2,10 +2,27 @@ function sol = F_par(x1,x2,y1,y2,d1,d2,d3) % 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,y1+d1,sqrt((x2-y2-d2)^2+d3^2))... - -(x2-y2-d2)*g0(0.5,x2,y2+d2,sqrt((x1-y1-d1)^2+d3^2))... - +1/3*((x1-y1-d1)^2+(x2-y2-d2)^2+d3^2)^(1.5); + sol = (x1-y1-d1)*(x2-y2-d2); + + if(sol~=0) + sol = sol * G00(-0.5,x1,x2,y1+d1,y2+d2,d3); + end + + if((x1-y1-d1)~=0) + sol = sol - (x1-y1-d1)*g0(0.5,x1,y1+d1,sqrt((x2-y2-d2)^2+d3^2)); + end + + if((x2-y2-d2)~=0) + sol = sol -(x2-y2-d2)*g0(0.5,x2,y2+d2,sqrt((x1-y1-d1)^2+d3^2)); + end + + sol = sol + 1/3*((x1-y1-d1)^2+(x2-y2-d2)^2+d3^2)^(1.5); + + +% sol = (x1-y1-d1)*(x2-y2-d2)*G00(-0.5,x1,x2,y1+d1,y2+d2,d3)... +% -(x1-y1-d1)*g0(0.5,x1,y1+d1,sqrt((x2-y2-d2)^2+d3^2))... +% -(x2-y2-d2)*g0(0.5,x2,y2+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)]))... diff --git a/src/G00.m b/src/G00.m index f181612..1112d9f 100644 --- a/src/G00.m +++ b/src/G00.m @@ -3,16 +3,21 @@ function sol = G00(p,y1,y2,x1,x2,l) % sol = nan; if(p==-1.5) if(l==0) - sol = -sqrt((y1-x1)^2+(y2-x2)^2)/((y1-x1)*(y2-x2)); + sol = -sqrt((y1-x1)^2+(y2-x2)^2)/((y1-x1)*(y2-x2)); else 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*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); + sol = 2*p*l^2*G00(p-1,y1,y2,x1,x2,l); + if((y1-x1)~=0) + sol = sol + (y1-x1)*g0(p,y2,x2,sqrt((y1-x1)^2+l^2)); + end + if((y2-x2)~=0) + sol = sol + (y2-x2)*g0(p,y1,x1,sqrt((y2-x2)^2+l^2)); + end + sol = sol / (2*p+2); else error 'no case for p defined' end diff --git a/src/test_solve.m b/src/test_solve.m index 5f2e0e3..99c5f23 100644 --- a/src/test_solve.m +++ b/src/test_solve.m @@ -3,7 +3,6 @@ N = size(elements); A = zeros(N); - % untere schranke s t obere schranke k l intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ... f(k1,k2,l1,l2)-f(k1,k2,l1,t2)-f(k1,k2,t1,l2)+f(k1,k2,t1,t2)... @@ -11,18 +10,12 @@ intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ... -f(s1,k2,l1,l2)+f(s1,k2,l1,t2)+f(s1,k2,t1,l2)-f(s1,l2,t1,t2)... +f(s1,s2,l1,l2)-f(s1,0,l1,t2)-f(s1,s2,t1,l2)+f(s1,s2,t1,t2); -% intF(@(s1,s2,k1,k2,t1,t2,l1,l2) F_par(s1,s2,k1,k2,t1,t2,l1,l2,d(1),d(2),d(3))... -% ,s1,s2,k1,k2,t1,t2,l1,l2); for j = 1:N - - for k = 1:N - - ej = coordinates(elements(j,:)',:); ek = coordinates(elements(k,:)',:); - d = ek(1,:) - ej(1,:) - + d = ek(1,:) - ej(1,:); + A(j,k) = 1/ (4*pi) *... intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))... ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2)); @@ -32,6 +25,4 @@ end A b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2)) - - -x = A\b; \ No newline at end of file +x = A\b \ No newline at end of file -- 2.47.3