ACCA-Imperial / SKPrime

A MATLAB implementation of the Schottky-Klein prime function.
GNU General Public License v3.0
10 stars 6 forks source link

greensC0 real constant not correctly evaluated #29

Closed rhodrin closed 8 years ago

rhodrin commented 8 years ago

When evaluating the function greensC0 at specific locations the 'undetermined real constant' from the Schwarz problem is not evaluated correctly. Solutions using this function differ from those in which G0 is formed using the Schottky Klein Prime function by a real constant (which for practical applications should not be the case).

Probably also applied to GJ (but have not checked) g0constant.zip

.

rhodrin commented 8 years ago

*applies to GJ

ehkropf commented 8 years ago

@rhodrin Do you know of a way to determine this constant a priori? Does this real value exist independent of the method used to find the Green's function?

rhodrin commented 8 years ago

Their should be, and I imagine it's related to the argument of alpha. I'll think about it this evening.

Another thing though (and I'm probably missing something here or have made a mistake), g0constant.zip but in the following m file why do the real parts of G0's values on the boundaries differ by a non-constant?

ehkropf commented 8 years ago

Try this instead to construct g0a:

w1 = skprime(alpha, D);
w1c = invParam(w1);
g0a = @(z) log(w1(z)./(abs(alpha).*w1c(z)))/(2i*pi);

(Note invParam gives the prime function with the parameter inverted wrt the unit circle, i.e. 1/conj(alpha).) I get the values differ by a real constant when I do this.

rhodrin commented 8 years ago

So I believe the constant can be evaluated by imposing a normalisation on \hat(G)_0 (not G_0). If \hat(G)_0 is defined as in eq. (5.9) of the paper, then as zeta -> alpha, \hat(G)_0=1/(2i_pi)_log(1/(abs(alpha)*\hat(w)(alpha,1/conj(alpha)))), where w(zeta,alpha)=(zeta-alpha)\hat(w)(zeta,alpha) is the SKprime function.

So that normalisation should work, but relies on knowing \hat(w)(alpha,1/conj(alpha)). I'll try and think of another that only relies on known properties of the Greens functions but am not certain it can be done.

ehkropf commented 8 years ago

Ok, cool, I'll have a look. I'm not too keen on needing the prime function to normalise G0(hat) however.

rhodrin commented 8 years ago

Ok got it! (I think).

Just evaluate \hat(G)_0 at alpha. The result should be purely complex. Therefore the real part that comes from evaluating g0.g0hatEval(alpha) is the real constant that needs to be factored out. (-0.0015 in the example below).

dv = [
    0.0-0.5i
    -0.1+0.5i];
qv = [
    0.1
    0.1];

D = skpDomain(dv, qv);

zp = boundaryPts(D, 10);
alpha=0.3*exp(0.2i*pi);

g0 = greensC0(alpha, D);

w1 = skprime(alpha, D);
w2 = invParam(w1);

Xh = @(z) w2.Xeval(z)./(z-1./conj(alpha)).^2;

g0ha = @(z) 1./(4i*pi).*log(1./(alpha*conj(alpha)*Xh(z)));

g0hn = @(z) g0.g0hatEval(z);

disp(g0ha(alpha))
disp(g0hn(alpha))
rhodrin commented 8 years ago

And here's some code illustrating it working.


%% G0 evaluated correctly.
clear

%%

dv = [
    0.0-0.5i
    -0.1+0.5i];
qv = [
    0.1
    0.1];

D = skpDomain(dv, qv);

zp = boundaryPts(D, 10);
alpha=0.3*exp(0.2i*pi);

%%

g0 = greensC0(alpha, D);

g0hn = @(z) g0.g0hatEval(z);

const=real(g0hn(alpha));

gf=g0(zp)-const;

% disp(gf)

%%

w1 = skprime(alpha, D);
w2 = invParam(w1);

g0a = @(z) 1./(2i*pi).*log(w1(z)./(abs(alpha).*w2(z)));

ga=g0a(zp);

% disp(ga)

disp(ga-gf)
ehkropf commented 8 years ago

Excellent. I will have a look.

ehkropf commented 8 years ago

Ugh. Breaks again with the domain

dv = [
  0.052448+0.36539i
  -0.27972-0.12762i
   0.48252-0.28147i];
qv = [
  0.15197
  0.17955
  0.20956];
rhodrin commented 8 years ago

What alpha are you seeing things break with? I've tried a few different ones and they seem to work.

rhodrin commented 8 years ago

I think for the above domain and alpha=0.3_exp(0.2i_pi) say, the 'g0 evaluation error' is kicking in and the problem doesn't lie with this fix?

ehkropf commented 8 years ago

Just saw the same thing and had that same thought. Testing it now.

ehkropf commented 8 years ago

It's the correction algorithm screwing things up, because g0.hat is no longer what you think it is. Checking if there is a way to compensate.

ehkropf commented 8 years ago

Oh, and try alpha = dv(1) + (qv(1) + dh)*exp(0.25i*pi); with dh=0.15 and dh=0.14. Any parameter within a distance 0.15 of a circle triggers the correction code.

ehkropf commented 8 years ago

It also bothers me this constant, when it works, isn't to machine accuracy.

rhodrin commented 8 years ago

Ok, will try now. I was just using this code and seeing different complex components for g0hat (for certain alpha) so knew something funny was going on.


clear

dv = [
  0.052448+0.36539i
  -0.27972-0.12762i
   0.48252-0.28147i];
qv = [
  0.15197
  0.17955
  0.20956];

D = skpDomain(dv, qv);

zp = boundaryPts(D, 10);
alpha=0.3*exp(0.2i*pi);

for j=1:3
    disp(abs(alpha-dv(j))-qv(j))
end

g0 = greensC0(alpha, D);

w1 = skprime(alpha, D);
w2 = invParam(w1);

Xh = @(z) w2.Xeval(z)./(z-1./conj(alpha)).^2;

g0ha = @(z) 1./(4i*pi).*log(1./(alpha*conj(alpha)*Xh(z)));

g0hn = @(z) g0.g0hatEval(z);

disp(g0ha(alpha))
disp(g0hn(alpha))
rhodrin commented 8 years ago

Yep, using alpha = dv(1) + (qv(1) + dh)_exp(0.25i_pi); with dh=0.15 in the above code the complex parts of g0ha and g0hn agree. When dh is changed to 0.14 they disagree.

rhodrin commented 8 years ago

This is the correction explained in appendix D right? I'll take a read over that.

ehkropf commented 8 years ago

Well that was easy. Try

const = real(g0.hat(alpha) + log(g0.singCorrFact(alpha))/(2i*pi));

(Yep, the correction factor is in the appendix.)

rhodrin commented 8 years ago

Cool.