Koch曲线的matlab程序实现

Koch曲线是一个典型的“数学怪物”，它的特点是处处连续但点点不可微。下面将讲解Koch曲线的形成原理及matlab程序实现（对空间上的任意两个点，画出它们的Koch曲线）

形成过程

对于二维空间上的任意一条线段，将该线段的两个端点称为首点和尾点，我们需要利用这首点和尾点来得到新的五个点。
首先，第一和第五个点即是首点和尾点，第二和第四个点是上述线段的两个三等分点，而第三个点是跟这两个三等分点组成等边三角形的三角形顶点。依次将五个点连接起来，便得到迭代一次的Koch曲线，然后对新得到的曲线上相邻的两个点做同样的操作。如此重复便可以得到Koch曲线，如下图所示。

对原始线段迭代一次之后得到的四折线段可以看成是原始线段分别做了四次变换之后得到的线段，一般的变换式为
Wi[xy]=ωi[xy]+ti(i=1,2,3,4)W_i \begin{bmatrix} x\\ y \\ \end{bmatrix}=\omega_i \begin{bmatrix} x\\ y \\ \end{bmatrix}+t_i(i=1,2,3,4) Wi[xy]=ωi[xy]+ti(i=1,2,3,4)
其中
ω1=[130013],t1=[00],\omega_1= \begin{bmatrix} \frac{1}{3} & 0 \\ 0 & \frac{1}{3} \\ \end{bmatrix}, t_1= \begin{bmatrix} 0 \\ 0 \\ \end{bmatrix}, ω1=[310031],t1=[00],
ω2=[16−363616],t2=[130],\omega_2= \begin{bmatrix} \frac{1}{6} & -\frac{\sqrt{3}}{6} \\ \frac{\sqrt{3}}{6} & \frac{1}{6} \\ \end{bmatrix}, t_2= \begin{bmatrix} \frac{1}{3} \\ 0 \\ \end{bmatrix}, ω2=[6163−6361],t2=[310],
ω3=[1636−3616],t3=[1236],\omega_3= \begin{bmatrix} \frac{1}{6} & \frac{\sqrt{3}}{6} \\ -\frac{\sqrt{3}}{6} & \frac{1}{6} \\ \end{bmatrix}, t_3= \begin{bmatrix} \frac{1}{2} \\ \frac{\sqrt{3}}{6} \\ \end{bmatrix}, ω3=[61−636361],t3=[2163],
ω4=[130013],t4=[230],\omega_4= \begin{bmatrix} \frac{1}{3} & 0 \\ 0 & \frac{1}{3} \\ \end{bmatrix}, t_4= \begin{bmatrix} \frac{2}{3}\\ 0 \\ \end{bmatrix}, ω4=[310031],t4=[320],

程序实现

定义Koch曲线函数

function Koch_line(time, A)
% time：迭代次数 A：初始两点坐标（列向量）
% time=10;A=[0,1; 0,0];
w1 = [1.0/3, 0; 0, 1.0/3]; t1 = [0, 0];
w2 = [1.0/6, -sqrt(3)/6; sqrt(3)/6, 1.0/6]; t2 = [1.0/3, 0];
w3 = [1.0/6, sqrt(3)/6; -sqrt(3)/6, 1.0/6]; t3 = [1.0/2, sqrt(3)/6];
w4 = [1.0/3, 0; 0, 1.0/3]; t4 = [2.0/3, 0];
for i=1:timetest_A = zeros(2, 1+4^i);[m,n] = size(A);for j=1:n-1 % 每次储存4个节点test_A(:, (j-1)*4+1) = A(:, j); % 首点test_A(:, (j-1)*4+2) = 1/3 * (A(:, j+1)-A(:, j)) + A(:, j); % 首点通过第二个变换得到test_A(:, (j-1)*4+3) = w2*(A(:, j+1)-A(:, j)) + A(:, j) + 1/3 * (A(:, j+1)-A(:, j)); % 尾点通过第二个变换得到test_A(:, (j-1)*4+4) = 2/3 * (A(:, j+1)-A(:, j)) + A(:, j); % 首点通过第四个变换得到endtest_A(:, end) = A(:, end);A = test_A;
end
figure,line(A(1, :), A(2, :))
axis([0 1 0 1])
end

演示

time=10;A=[0,1; 0,0];
Koch_line(time, A)

time=10;A=[0,1; 0,1];
Koch_line(time, A)

注意

1、点的计算

已知首点和尾点，现在需要得到新的五个点。
首先，第一和第五个点就是原先的首尾点，其次，第二个点可以由尾点通过变换W1W_1W1得到，也可以由首点通过变换W2W_2W2得到，第三个点可以由尾点通过变换W2W_2W2得到，也可以由首点通过变换W3W_3W3得到，第四个点可以由尾点通过变换W3W_3W3得到，也可以由首点通过变换W4W_4W4得到。
所以我在写程序的时候第一和第五个直接由首尾点得到，第二个点是由首点通过变换W2W_2W2得到，第三个点是由尾点通过变换W2W_2W2得到，第四个点是由首点通过变换W4W_4W4得到，这样就能减少使用变换W3W_3W3。

2、是否绕原点旋转

通过变换Wi,(i=1,2,3,4)W_i,(i=1,2,3,4)Wi,(i=1,2,3,4)，可以用首点和尾点的坐标得到新的五个点的坐标，但同时需要注意到，这种变换只针对原来的两个点有一个是原点的情况下，其中像变换W2W_2W2和W3W_3W3，它们是将线段绕原点进行旋转，再加上一个平移量。而对于空间上的任意两个点进行变换，如果仍然采用上述的变换矩阵，将不会得到希望的结果。
这时就需要对首点和尾点进行平移，将它们的旋转中心平移到原点，再进行旋转变换，之后再平移回去。
举个例子能更明白一些，比如对首尾点相连的线段lll进行变换W2W_2W2，首先是将lll绕首点逆时针旋转60度，然后将lll缩小为原来长度的三分之一，最后沿尾点的方向平移原来lll长度的三分之一距离。如果首点不是原点，则进行上述的变换并不能得到所要的结果，这时需要先将首尾点的坐标值减去首点的坐标值，在进行上述变换，最后再把首点的坐标值加上，如此操作才能画出正确的Koch曲线。

----end----

参考资料
李水根，吴纪桃《分形与小波》