Instead of the proposed integration over the area under the Koch curve, it's easier
to employ Green's theorem and integrate over the x-axis and then backwards over the
Koch curve (i.e. in counterclockwise order):
$$
\iint \left[ \frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y} \right] \,dx\,dy =
\oint \left[ P\,dx + Q\,dy \right]
$$
There are several possible choices for $P(x,y)$ and $Q(x,y)$ but this one will be ours:
$$
P(x,y) = -\frac{e^{a \cdot x + b \cdot y}}{b} \quad ; \quad Q(x,y) = 0
$$
Giving:
$$
\iint_K e^{a \cdot x + b \cdot y} \,dx\,dy = - \frac{1}{b} \oint_K e^{a \cdot x + b \cdot y}\,dx =
- \frac{1}{b}\int_0^1 e^{a \cdot x + 0}\,dx - \frac{1}{b} \int_K e^{a \cdot x + b \cdot y}\,dx
$$
Where the last integral is over the Koch curve from the right $(1,0)$ to the left $(0,0)$. So:
$$
\iint_K e^{a \cdot x + b \cdot y} \,dx\,dy =
\left(1-e^a\right)/(ab) + \frac{1}{b}\int_K e^{a \cdot x + b \cdot y}\,dx
$$
If we reverse the integration over the Koch curve from the left $(0,0)$ to the right $(1,0)$.
The method as proposed in the question is thus replaced by the following:
- Calculate the integral for an arbitrary line segment with
vertex coordinates $(x_1,y_1),(x_2,y_2)$
- Find the vertices of all line segments in a finite iterand of the Koch curve
- Find the integral for this finite interand, which is the sum of known integrals
over line segments
- Let the number of line segments grow to infinity; determine the limit eventually
(and/or approximately)
First we must calculate the exponent for an arbitrary line segment $S$.
Let
$\,x = x_1 + (x_2-x_1)t\,$ and $\,y = y_1 + (y_2-y_1)t$ , then:
$$
\int_{S} e^{a \cdot x+ b \cdot y}\, dx =
\int_0^1 e^{a \left[x_1 + (x_2-x_1)t \right]
+ b \left[y_1 + (y_2-y_1)t \right]}\;(x_2-x_1)\,dt =\\
(x_2-x_1)\frac{e^{a x_2 + b y_2}-e^{a x_1 + b y_1}}{(a x_2 + b y_2)-(a x_1 + b y_1)}
$$
The following (Pascal) program defines finite iterands of a Koch curve,
recursively. Arbitrary values of $a$ and $b$ can be substituted in the
test procedure eventually. Several tests are built in to ensure that
the program works correctly. The variable
sum is used to sum up for
the requested
Zach integral, but another variable
som is used to sum up for
the
Area under the Koch curve, which is exactly known for each iterand.
The outcome of the integral for one triangle, as mentioned in the question,
can be used for verifying the current integration procedure for the first
interand
Koch(1,a,b,1).
program Zachter;
type
vertex = record
x,y : double;
end;
segment = array[1..2] of vertex;
triangle = array[1..3] of vertex;
var
phi,p,q,sum,som : double;
function opp(S : segment; a,b : double) : double;
var
x_1,x_2, y_1,y_2, f : double;
begin
x_1 := S1.x; x_2 := S[2].x;
y_1 := S1.y; y_2 := S[2].y;
f := (y_1+y_2)/2;
opp := f*(x_2-x_1);
end;
function lijn(S : segment; a,b : double) : double;
var
x_1,x_2, y_1,y_2, f : double;
begin
x_1 := S1.x; x_2 := S[2].x;
y_1 := S1.y; y_2 := S[2].y;
f := (exp(a*x_2+b*y_2)-exp(a*x_1+b*y_1))
/ ((a*x_2+b*y_2)-(a*x_1+b*y_1));
lijn := f*(x_2-x_1)/b;
end;
procedure voorts(len : double; a,b : double);
var
S : segment;
begin
S1.x := p; S1.y := q;
p := p+cos(phi)*len;
q := q+sin(phi)*len;
S[2].x := p; S[2].y := q;
sum := sum + lijn(S,a,b);
som := som + opp(S,a,b);
end;
procedure turnleft(angle : integer);
begin
phi := phi+angle*Pi/180;
end;
procedure turnright(angle : integer);
begin
phi := phi-angle*Pi/180;
end;
procedure Koch(len,a,b: Double; count: Integer);
begin
If count = 0 then
voorts(len,a,b)
else begin
Koch(len/3,a,b, count-1);
turnleft(60);
Koch(len/3,a,b, count-1);
turnright(120);
Koch(len/3,a,b, count-1);
turnleft(60);
Koch(len/3,a,b, count-1);
end;
end;
function kernel(D : triangle; a,b : double) : double;
var
x_1,x_2,x_3, y_1,y_2,y_3, f : double;
begin
x_1 := D1.x; x_2 := D[2].x; x_3 := D[3].x;
y_1 := D1.y; y_2 := D[2].y; y_3 := D[3].y;
f := -(exp(a*x_2+b*y_2)*a*x_1+exp(a*x_2+b*y_2)*b*y_1-exp(a*x_2+b*y_2)
*a*x_3-exp(a*x_2+b*y_2)*b*y_3+exp(a*x_3+b*y_3)*a*x_2-exp(a*x_3+b*y_3)
*a*x_1+exp(a*x_3+b*y_3)*b*y_2-exp(a*x_3+b*y_3)*b*y_1-exp(a*x_1+b*y_1)
*a*x_2+exp(a*x_1+b*y_1)*a*x_3-exp(a*x_1+b*y_1)*b*y_2+exp(a*x_1+b*y_1)
*b*y_3)/(-a*x_3+a*x_1-b*y_3+b*y_1)/(a*x_2-a*x_3+b*y_2-b*y_3)
/(-a*x_2+a*x_1-b*y_2+b*y_1);
kernel := f*((x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1));
end;
procedure test(a,b : double);
var
k : integer;
f,t : double;
D : triangle;
begin
f := 4/9; t := 1;
for k := 0 to 12 do
begin
p := 0; q := 0; phi := 0;
som := 0;
sum := (1-exp(a))/(a*b);
Koch(1,a,b,k);
{ http://nl.wikipedia.org/wiki/Koch-kromme }
Writeln('Area =',som,' =',(1-t)*9/5*sqrt(3)/36);
if k=1 then
begin
Write('Zach =',sum,'= ');
D1.x := 1/3; D1.y := 0;
D[2].x := 2/3; D[2].y := 0;
D[3].x := 1/2; D[3].y := sqrt(3)/6;
Writeln(kernel(D,a,b));
end else Writeln('Zach =',sum);
t := t*f;
end;
end;
begin
test(2,3);
end.
Output for $a=2$ and $b=3$ , thirteen interations:
Area = 0.00000000000000E+0000 = 0.00000000000000E+0000
Zach = 0.00000000000000E+0000
Area = 4.81125224324688E-0002 = 4.81125224324688E-0002
Zach = 1.79819984312336E-0001= 1.79819984312336E-0001
Area = 6.94958657357883E-0002 = 6.94958657357883E-0002
Zach = 2.67676463467801E-0001
Area = 7.89995738705970E-0002 = 7.89995738705969E-0002
Zach = 3.07390573739620E-0001
Area = 8.32234441527343E-0002 = 8.32234441527341E-0002
Zach = 3.25080441215073E-0001
Area = 8.51007198336841E-0002 = 8.51007198336840E-0002
Zach = 3.32944604962222E-0001
Area = 8.59350645807717E-0002 = 8.59350645807728E-0002
Zach = 3.36439888372271E-0001
Area = 8.63058844683464E-0002 = 8.63058844683678E-0002
Zach = 3.37993352588789E-0001
Area = 8.64706933071906E-0002 = 8.64706933072990E-0002
Zach = 3.38683781372148E-0001
Area = 8.65439416795354E-0002 = 8.65439416801572E-0002
Zach = 3.38990638617374E-0001
Area = 8.65764965093825E-0002 = 8.65764965125387E-0002
Zach = 3.39127019600374E-0001
Area = 8.65909653108997E-0002 = 8.65909653269305E-0002
Zach = 3.39187633284637E-0001
Area = 8.65973958352309E-0002 = 8.65973959111046E-0002
Zach = 3.39214572284906E-0001
The outcome is expected to converge, but not very fast,
and there is quite some loss of precision, as is seen
with the Area outputs.
