课后作业3

Q1. 三个分形图形生成的程序解释 每行代码有什么用

newton.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
function seta=f_arg(x,y)
% 定义一个名为 f_arg 的函数，输入为 x 和 y，输出为 seta (theta)

if x==0
% 情况1：如果 x=0，点在 y 轴上
    if y>=0
    % 如果 y >= 0 (点在 y 轴正半轴或原点)
        seta=pi/2; % 角度为 90 度 (pi/2)
    else
    % 如果 y < 0 (点在 y 轴负半轴)
        seta=-pi/2; % 角度为 -90 度 (-pi/2)
    end
elseif x>0
% 情况2：如果 x>0，点在右半平面 (第一、四象限)
    seta=atan(y/x); % 直接使用反正切函数 atan 即可计算出正确的角度
elseif x<0 && y>=0
% 情况3：如果 x<0 且 y>=0，点在左上方的第二象限
    seta=pi+atan(y/x); % atan(y/x) 的结果在 (-pi/2, 0] 之间，需要加上 pi 才能映射到正确的角度
elseif x<0 && y<0
% 情况4：如果 x<0 且 y<0，点在左下方的第三象限
    seta=-pi+atan(y/x); % atan(y/x) 的结果在 (0, pi/2) 之间，需要减去 pi 才能映射到正确的角度
else %修正2：原始分支条件仍不完备，bug时x,y=NaN
% 情况5 (特殊情况)：注释表明这是为了处理 x 或 y 是 NaN (Not a Number) 的情况
    seta=pi/2; % 当出现 bug 时，按原点附近的特殊情况处理，给一个默认值
end
end

function f=f_iterat(x,y,n,m,delta)
% 定义函数 f_iterat，输入初始点(x,y)，方程次数n，最大迭代次数m，收敛判据delta
% 输出 f 是收敛时的迭代次数，如果不收敛则为 0

f=0; % 初始化返回值为 0，表示默认不收敛
for k=1:m
% 开始一个循环，最多迭代 m 次，k 是当前迭代的次数
    
    % --- 将笛卡尔坐标 (x,y) 转换为极坐标 (r, seta) ---
    r=sqrt(x^2+y^2); % 计算当前点的模长 r (到原点的距离)。注意：原文的 xx+yy 是笔误，应为 x^2+y^2。
    if r==0 %修正1：为避免出现inf值，强制给个极小值
    % 如果点在原点，r=0 会导致后续计算出现除以0的错误
        r=1e-9; % 将 r 设为一个非常小的正数来避免错误
    end
    seta=f_arg(x,y); % 调用 f_arg 函数计算当前点的辐角
    
    % --- 牛顿法迭代公式 z_new = z - f(z)/f'(z) 的极坐标形式 ---
    % 对于方程 f(z) = z^n - 1，其牛顿迭代公式为 z_k+1 = ((n-1)z_k^n + 1) / (n * z_k^(n-1))
    % 将 z = r * (cos(seta) + i*sin(seta)) 代入并化简，就得到下面的 x2 和 y2 的表达式
    x2=((n-1)*r*cos(seta)+r^(1-n)*cos((1-n)*seta))/n; % 计算下一次迭代点的 x 坐标。注意：原文 rcos 是笔误，应为 r*cos。
    y2=((n-1)*r*sin(seta)+r^(1-n)*sin((1-n)*seta))/n; % 计算下一次迭代点的 y 坐标。注意：原文 rsin 是笔误，应为 r*sin。
    
    % --- 判断是否收敛 ---
    dis=sqrt((x2-x)^2+(y2-y)^2); % 计算新点 (x2,y2) 和旧点 (x,y) 之间的距离
    if dis<delta
    % 如果这个距离小于预设的阈值 delta，我们认为迭代已经收敛
        f=k; % 将当前的迭代次数 k 赋给返回值 f
        break % 跳出 for 循环，因为已经找到了结果
    end
    
    % --- 更新迭代点 ---
    x=x2; % 将新点的坐标赋给旧点，准备下一次循环
    y=y2;
end
end

function newton(n)%牛顿迭代法产生分形图片
% 定义主函数 newton，输入参数 n 代表方程 z^n - 1 = 0 的次数

clc % 清空 MATLAB 的命令行窗口
m=16; % 设置最大迭代次数为 16
a1=-3; % 设置绘图区域的 x 轴左边界
a2=3;  % 设置绘图区域的 x 轴右边界
b1=-2; % 设置绘图区域的 y 轴下边界
b2=2;  % 设置绘图区域的 y 轴上边界
delta=0.001; % 设置收敛判断的阈值为 0.001
hold on % 告诉 MATLAB 接下来的所有 plot 命令都在同一张图上绘制，而不是新建或覆盖
for x=a1:0.02:a2
% 外层循环：遍历 x 坐标。从 a1(-3) 到 a2(3)，步长为 0.02
    for y=b1:0.02:b2
    % 内层循环：遍历 y 坐标。从 b1(-2) 到 b2(2)，步长为 0.02
    % 这两个循环组合起来，就构成了一个覆盖整个绘图区域的网格
    
        f=f_iterat(x,y,n,m,delta); % 对网格上的每一个点 (x,y) 调用迭代函数，获取收敛速度 f
        
        if f>0
        % 如果 f > 0，说明这个点在 m 次迭代内收敛了
            switch rem(f,7)
            % 使用 switch 语句根据收敛速度 f 来选择颜色
            % rem(f,7) 是计算 f 除以 7 的余数，结果范围是 0 到 6
            % 这样可以将收敛速度相近的点映射到同一种颜色，形成分形的彩色条带
                case 0
                    plot(x,y,'w.'); % 余数为0，绘制白色(white)的点
                case 1
                    plot(x,y,'b.'); % 余数为1，绘制蓝色(blue)的点
                case 2
                    plot(x,y,'r.'); % 余数为2，绘制红色(red)的点
                case 3
                    plot(x,y,'y.'); % 余数为3，绘制黄色(yellow)的点
                case 4
                    plot(x,y,'g.') ; % 余数为4，绘制绿色(green)的点
                case 5
                    plot(x,y,'c.') ; % 余数为5，绘制青色(cyan)的点
                case 6
                    plot(x,y,'m.'); % 余数为6，绘制品红色(magenta)的点
            end
        end
    end
end
hold off % 绘图结束，关闭 hold on 状态，后续的 plot 命令会新建图形
end

这是一个用于生成牛顿分形（Newton Fractal） 图像的 MATLAB 程序。牛顿分形是通过牛顿迭代法求解复数域上的方程 z^n - 1 = 0 时，根据初始点的不同，迭代收敛到不同根的速度来着色的。

整个程序由三个函数组成：

newton(n): 主函数，负责设置绘图区域、遍历像素点并调用迭代函数，最后根据结果绘图。
f_iterat(…): 核心迭代函数，对给定的初始点 (x, y) 进行牛顿迭代，并返回收敛所需的迭代次数。
f_arg(x,y): 辅助函数，用来计算一个点（复数）的辐角（angle），即 atan2(y,x) 的功能。

juliaSet.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
phi = inline('z^2 - 1.25');   % 定义内联函数phi(z) = z² - 1.25，用于寻找其不动点（满足phi(z)=z的z值）
fixpt1 = (1 + sqrt(6))/2;     % 计算函数phi的第一个不动点：解方程z² - 1.25 = z（即z² - z - 1.25 = 0），根为(1+√6)/2
fixpt2 = (1 - sqrt(6))/2;     % 计算函数phi的第二个不动点：方程的另一个根(1-√6)/2

colormap([1 0 0; 1 1 1]);     % 设置颜色映射：索引1对应红色（RGB[1,0,0]），索引2对应白色（RGB[1,1,1]）
                              % 后续用1标记"有界轨道点"（红色），2标记"无界轨道点"（白色）
M = 2*ones(141,361);          % 初始化141×361的矩阵M，所有元素为2（默认标记为白色，即初始假设所有点轨道无界）

for j=1:141,                  % 外循环：遍历虚部范围（j控制虚部）
  y = -.7 + (j-1)*.01;        % 计算当前点的虚部y：范围从-0.7到0.7（步长0.01，共141个点：(0.7 - (-0.7))/0.01 + 1 = 141）
  for i=1:361,                % 内循环：遍历实部范围（i控制实部）
    x = -1.8 + (i-1)*.01;     % 计算当前点的实部x：范围从-1.8到1.8（步长0.01，共361个点：(1.8 - (-1.8))/0.01 + 1 = 361）
    z = x + 1i*y;             % 构造复数初始点z（实部x，虚部y），1i是MATLAB中表示虚数单位√(-1)的符号
    zk = z;                   % 初始化迭代变量zk为初始点z
    iflag1 = 0;               % 标志1：计数连续迭代中接近第一个不动点fixpt1的次数（初始为0）
    iflag2 = 0;               % 标志2：计数连续迭代中接近第二个不动点fixpt2的次数（初始为0）
    kount = 0;                % 迭代计数器：记录总迭代次数（初始为0）

    % 循环条件：迭代次数<100、当前迭代值zk的模<2（避免发散到无穷）、连续接近fixpt1的次数<5、连续接近fixpt2的次数<5
    while kount < 100 & abs(zk) < 2 & iflag1 < 5 & iflag2 < 5,
      kount = kount+1;        % 迭代次数加1
      zk = phi(zk);           % 执行不动点迭代：zk更新为phi(zk) = zk² - 1.25
      err1 = abs(zk-fixpt1);  % 计算当前zk与第一个不动点的距离（误差）
      if err1 < 1.e-6,        % 若距离小于1e-6（认为接近fixpt1）
         iflag1 = iflag1 + 1; % 连续接近次数加1
      else
         iflag1 = 0;          % 否则重置为0（中断连续计数）
      end;
      err2 = abs(zk-fixpt2);  % 计算当前zk与第二个不动点的距离（误差）
      if err2 < 1.e-6,        % 若距离小于1e-6（认为接近fixpt2）
        iflag2 = iflag2 + 1;  % 连续接近次数加1
      else
        iflag2 = 0;           % 否则重置为0
      end;
    end;
    % 若连续5次接近某个不动点（收敛），或迭代满100次（可能有界但未收敛），则标记为红色（1）
    if iflag1 >= 5 | iflag2 >= 5 | kount >= 100,   
      M(j,i) = 1;                                  
    end;
  end;
end;

image([-1.8 1.8],[-.7 .7],M),  % 绘制图像：x轴范围[-1.8,1.8]，y轴范围[-0.7,0.7]，用矩阵M的数值（1/2）索引颜色映射
axis xy                        % 设置坐标轴方向：y轴向上为正（默认image函数y轴向下，需纠正）

代码核心逻辑是：对复平面上特定区域（实部 [-1.8,1.8]、虚部 [-0.7,0.7]）内的大量初始点，通过迭代函数phi(z)=z²-1.25生成轨道，判断轨道是否有界（收敛到不动点或迭代 100 次仍未发散）。有界点标记为红色，无界点（快速发散）标记为白色，最终绘制出该迭代过程的分形边界图案

‍

yoda.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
%% Load the model/tessellation information  % 代码段标题：加载模型/网格信息
load yodapose_low  % 加载名为yodapose_low的.mat文件，包含模型的顶点(V)、面(F3/F4)等数据

%% Create initial plot  % 代码段标题：创建初始图形
Vt = V;  % 将原始顶点数据V复制到Vt（Vt用于后续动画中的顶点变换，保留原始V不变）
clf  % 清除当前图形窗口内容
% 绘制3D模型：使用顶点Vt和面F3，设置面颜色为浅绿色[0.76, 0.87, 0.78]
patch('Vertices',Vt,'Faces',F3,'FaceColor',[.76 .87 .78]);
% 绘制模型的另一部分面F4，使用相同颜色（F3和F4可能是模型不同区域的面集合）
patch('Vertices',Vt,'Faces',F4,'FaceColor',[.76 .87 .78]);
% 设置坐标轴：紧凑显示(tight)、各轴比例相等(equal)、保持3D视角(vis3d，避免旋转时视角突变)
axis tight equal vis3d
drawnow  % 强制更新图形，立即显示初始模型

%% Create translation matrix  % 代码段标题：创建平移矩阵（后续平移动画的参数，当前被注释的动画使用）
slides = 24;  % 动画的总帧数（平移动画计划分24帧完成）

% Create the translation matrix  % 创建平移矩阵
yMinValue = min(V(:,2,:));  % 提取模型所有顶点的y坐标（假设V的第二列是y轴），计算最小值
axisValues = axis;  % 获取当前坐标轴范围，返回格式为[xmin xmax ymin ymax zmin zmax]
yAxesMax = axisValues(4);  % 从坐标轴范围中提取y轴最大值（第4个元素）
% 计算每次平移的步长：总位移为(y轴最大值 - 顶点最小y值)，平均分配到24帧，确保最后一帧模型移出y轴范围
shift = (yAxesMax - yMinValue)/slides;  
[n,m] = size(V);  % 获取顶点矩阵V的尺寸：n是顶点数量，m=3（x、y、z三个维度）
% 创建平移矩阵T：n行3列，x和z方向平移量为0，y方向平移量为shift（每帧沿y轴移动shift）
T = [zeros(n,1),shift*ones(n,1),zeros(n,1)];  

% %% Animate translation  % 代码段标题：平移动画（当前被注释，未执行）
% 
% for i=1:slides  % 循环24帧，执行平移动画
%     Vt = Vt + T;  % 顶点沿y轴平移shift（累计位移随帧数增加）
%     cla  % 清除当前坐标轴内容（避免重复绘制）
%     % 重新绘制平移后的模型（面F3和F4）
%     patch('Vertices',Vt,'Faces',F3,'FaceColor',[.76 .87 .78]);
%     patch('Vertices',Vt,'Faces',F4,'FaceColor',[.76 .87 .78]);
%     axis(axisValues)  % 保持坐标轴范围不变（与初始范围一致）
%     drawnow  % 强制更新图形，显示当前帧
%     hold on  % 保持当前绘图（此处作用不大，因cla已清除）
%     pause(0.2)  % 暂停0.2秒，控制动画速度
% end

theta=0:pi/24:pi/2;  % 生成旋转角度序列：从0到π/2（0到90度），步长π/24，共13个角度（含首尾）
n=length(theta);  % 获取角度序列的长度（即旋转动画的总帧数）
for i=1:n  % 循环每一个角度，执行旋转动画
    % 定义绕y轴旋转的旋转矩阵R（右手坐标系）：
    % 第一行：[cosθ, 0, -sinθ]（x和z坐标受旋转影响）
    % 第二行：[0, 1, 0]（y坐标不变，因绕y轴旋转）
    % 第三行：[sinθ, 0, cosθ]
    R=[cos(theta(i)) 0 -sin(theta(i));0 1 0;sin(theta(i)) 0 cos(theta(i))];
    Vt = Vt*R;  % 顶点矩阵与旋转矩阵相乘，得到旋转后的顶点坐标（Vt每行是一个顶点的[x,y,z]，右乘旋转矩阵实现旋转）
    cla  % 清除当前坐标轴内容
    % 重新绘制旋转后的模型（面F3和F4）
    patch('Vertices',Vt,'Faces',F3,'FaceColor',[.76 .87 .78]);
    patch('Vertices',Vt,'Faces',F4,'FaceColor',[.76 .87 .78]);
    axis(axisValues)  % 保持坐标轴范围与初始一致
    drawnow  % 强制更新图形，显示当前旋转帧
    hold on  % 保持当前绘图（作用同上）
    pause(0.2)  % 暂停0.2秒，控制动画速度
end

代码核心功能是：加载一个 3D 模型（包含顶点 V 和两个面集合 F3、F4），先显示模型初始状态，然后通过循环实现模型的旋转动画（绕 y 轴从 0° 旋转到 90°，分 13 帧，每帧暂停 0.2 秒）。

‍

Q2. 根据三种方法创造自己的分形图形

newton

只需要取一个迭代次数n即可，这里取n=10n=10

代码：

1
2
3
clear
n = 10;
newton(n);

结果：

juliaSet

要想创造新的分形图形只需要修改迭代公式即可，这里将迭代公式定为：

z2+0.285+0.01iz^2+0.285+0.01i

则通过下式求其不动点

z2−x+0.285+0.01i=0z^2-x+0.285+0.01i=0

最终代码为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
phi = inline('z^2 + 0.285 + 0.01i');   % Define the function whose fixed points we seek.
fixpt1 = (1 + 4*(0.285 + 0.01i))/2;     % These are the fixed points.
fixpt2 = (1 - 4*(0.285 + 0.01i))/2;

colormap([1 0 0; 1 1 1]);     % Points numbered 1 (inside) will be colored red;
                              %   those numbered 2 (outside) will be colored white.
M = 2*ones(141,361);          % Initialize array of point colors to 2 (white).

for j=1:141,                  % Try initial values with imaginary parts between
  y = -.7 + (j-1)*.01;        %   -0.7 and 0.7
  for i=1:361,                % and with real parts between
    x = -1.8 + (i-1)*.01;     %   -1.8 and 1.8.
    z = x + 1i*y;             % 1i is the MATLAB symbol for sqrt(-1).
    zk = z;
    iflag1 = 0;               % iflag1 and iflag2 count the number of iterations
    iflag2 = 0;               %   when a root is within 1.e-6 of a fixed point;
    kount = 0;                % kount is the total number of iterations.

    while kount < 100 & abs(zk) < 2 & iflag1 < 5 & iflag2 < 5,
      kount = kount+1;
      zk = phi(zk);           % This is the fixed point iteration.
      err1 = abs(zk-fixpt1);  % Test for convergence to fixpt1.
      if err1 < 1.e-6,
         iflag1 = iflag1 + 1;
      else
         iflag1 = 0;
      end;
      err2 = abs(zk-fixpt2);  % Test for convergence to fixpt2.
      if err2 < 1.e-6,
        iflag2 = iflag2 + 1;
      else
        iflag2 = 0;
      end;
    end;
    if iflag1 >= 5 | iflag2 >= 5 | kount >= 100,   % If orbit is bounded, set this
      M(j,i) = 1;                                  %   point color to 1 (red).
    end;
  end;
end;

image([-1.8 1.8],[-.7 .7],M),  % This plots the results.
axis xy                        % If you don't do this, vertical axis is inverted.

结果：

yoda

这里通过python将奶龙的obj格式3D文件转换为mat文件

python程序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
import numpy as np
from scipy.io import savemat
import pywavefront

def obj_to_mat(obj_file_path, mat_file_path):
    """
    将OBJ文件转换为包含V、F3、F4的MAT文件（匹配示例格式）
    """
    # 读取OBJ文件
    scene = pywavefront.Wavefront(
        obj_file_path,
        collect_faces=True,
        parse=True
    )
    
    # 初始化数据（匹配目标格式）
    V = []         # 顶点坐标（N×3）
    F3 = []        # 三角形面（M×3，顶点索引）
    F4 = []        # 四边形面（K×4，顶点索引）
    
    # 提取所有顶点（OBJ中顶点是全局的，存储在scene.vertices）
    # 注意：pywavefront中顶点直接存储在scene.vertices（而非mesh.vertices）
    if hasattr(scene, 'vertices') and scene.vertices:
        V = np.array(scene.vertices)  # 转换为N×3矩阵
    
    # 遍历所有mesh提取面（区分三角形和四边形）
    for mesh in scene.mesh_list:
        if hasattr(mesh, 'faces') and mesh.faces:
            for face in mesh.faces:
                # OBJ面索引是1-based，转换为0-based（若需要1-based可删除-1）
                face_indices = np.array(face) + 1
                # 根据面的顶点数量分类
                if len(face_indices) == 3:
                    F3.append(face_indices)  # 三角形面
                elif len(face_indices) == 4:
                    F4.append(face_indices)  # 四边形面
                # （可选：处理更多边的多边形，这里忽略）
    
    # 转换为numpy矩阵（空则不保存）
    data = {}
    if V.size > 0:
        data['V'] = V  # 顶点坐标
    
    if F3:
        data['F3'] = np.array(F3)  # 三角形面
    
    if F4:
        data['F4'] = np.array(F4)  # 四边形面
    
    # 保存为MAT文件
    savemat(mat_file_path, data)
    print(f"成功转换：{obj_file_path} -> {mat_file_path}")
    print(f"包含的数据：{list(data.keys())}（V:顶点, F3:三角形面, F4:四边形面）")

# 使用示例
if __name__ == "__main__":
    obj_file = "nailong1.obj"    # 替换为你的OBJ文件路径
    mat_file = "nailong.mat"   # 输出MAT文件路径
    obj_to_mat(obj_file, mat_file)

通过修改变换矩阵，使得3D模型绕x轴旋转，由于这个模型的面数比较多，需要将面边缘的线画的更细，并加上一定的透明度，来获得更好的显示效果。

最终代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
clear
shg
%% 1. 加载模型数据
load nailong

%% 2. 初始化图形窗口
Vt = V;
clf;
hold on;
axis equal vis3d;  % 保持坐标轴比例，开启3D视角
patch('Vertices',V,'Faces',F3,'FaceColor',[.76 .87 .78],'EdgeColor','k','EdgeAlpha',0.1,'LineWidth',0.1);
drawnow;

%% 3. 计算平移参数（不变）
slides = 24;
yMinValue = min(V(:,2));
shift = (max(V(:,2)) - yMinValue + 0.1*max(V(:,2))) / slides;  % 适当放大平移量，避免边缘紧贴
n = size(V,1);
T = [zeros(n,1), shift*ones(n,1), zeros(n,1)];

% %% 4. 执行平移动画（取消固定坐标轴）
% for i = 1:slides
%     Vt = Vt + T;
%     cla;
%     patch('Vertices',Vt_rot,'Faces',F3,'FaceColor',[.76 .87 .78],'EdgeColor','k','EdgeAlpha',0.1,'LineWidth',0.1);
%     axis tight;  % 自动调整坐标轴，完整显示当前模型
%     axis equal;  % 保持比例不变
%     drawnow;
%     pause(0.2);
% end

%% 5. 执行上下旋转
theta = 0:pi/6:2*pi;
n_theta = length(theta);
for i = 1:n_theta
    % 绕x轴旋转矩阵
    R = [1          0           0;
        0    cos(theta(i))  -sin(theta(i));
        0    sin(theta(i))   cos(theta(i))];
    Vt_rot = Vt * R;
    cla;
    patch('Vertices',Vt_rot,'Faces',F3,'FaceColor',[.76 .87 .78],'EdgeColor','k','EdgeAlpha',0.1,'LineWidth',0.1);
    axis tight;  % 自动适配当前模型范围
    axis equal;  % 保持比例
    drawnow;
    pause(0.2);
end

hold off;

结果：

Q3. 牛顿迭代法产生分形图形时给大家的程序是用极坐标的形式来求复数根，请把牛顿迭代求复数跟的极坐标形式推导出来

这个分形图形是通过求解复数方程 zn−1=0z^n - 1 = 0 的根而产生的。这个方程在复数平面上有 n 个根，它们均匀分布在单位圆上。

牛顿迭代法的公式为：

zk+1=zk−f(zk)f′(zk)z_{k+1} = z_k - \frac{f(z_k)}{f'(z_k)}

其中，zkz_k 是第 k 次迭代的近似解，zk+1z_{k+1} 是下一次迭代的解。

我们要求解的方程是 zn−1=0z^n - 1 = 0。

所以，我们的目标函数是：

f(z)=zn−1f(z) = z^n - 1

对它求导，得到：

f′(z)=nzn−1f'(z) = n z^{n-1}

将 f(z)f(z) 和 f′(z)f'(z) 代入基本公式：

zk+1=zk−zkn−1nzkn−1z_{k+1} = z_k - \frac{z_k^n - 1}{n z_k^{n-1}}

⇒zk+1=nzkn−zkn+1nzkn−1\Rightarrow z_{k+1}=\frac{nz_{k}^{n}-z_{k}^{n}+1}{nz_{k}^{n-1}}

⇒zk+1=(n−1)zkn+1nzkn−1\Rightarrow z_{k+1}=\frac{(n-1)z_{k}^{n}+1}{nz_{k}^{n-1}}

⇒zk+1=(n−1)zknnzkn−1+1nzkn−1\Rightarrow z_{k+1}=\frac{(n-1)z_{k}^{n}}{nz_{k}^{n-1}}+\frac{1}{nz_{k}^{n-1}}

⇒zk+1=n−1nzk+1nzk1−n\Rightarrow z_{k+1}=\frac{n-1}{n}z_k+\frac{1}{n}z_{k}^{1-n}

设当前点 zk=x+iyz_k = x + iy。其极坐标形式为：

zk=r(cos⁡θ+isin⁡θ)z_k = r(\cos\theta + i\sin\theta)

其中，r=x2+y2r = \sqrt{x^2 + y^2} 是模长，θ=arg⁡(zk)\theta = \arg(z_k) 是辐角。

根据棣莫弗定理（De Moivre’s formula），复数的幂在极坐标下计算非常方便：

zkm=[r(cos⁡θ+isin⁡θ)]m=rm(cos⁡(mθ)+isin⁡(mθ))z_k^m = [r(\cos\theta + i\sin\theta)]^m = r^m(\cos(m\theta) + i\sin(m\theta))

我们将这个定理应用到迭代公式 zk+1=n−1nzk+1nzk1−nz_{k+1} = \frac{n-1}{n} z_k + \frac{1}{n} z_k^{1-n} 的第二项：

zk1−n=r1−n(cos⁡((1−n)θ)+isin⁡((1−n)θ))z_k^{1-n} = r^{1-n}(\cos((1-n)\theta) + i\sin((1-n)\theta))

现在，把 zkz_k 和 zk1−nz_k^{1-n} 的极坐标形式代回到迭代公式中：

zk+1=n−1n[r(cos⁡θ+isin⁡θ)]+1n[r1−n(cos⁡((1−n)θ)+isin⁡((1−n)θ))]z_{k+1} = \frac{n-1}{n} [r(\cos\theta + i\sin\theta)] + \frac{1}{n} [r^{1-n}(\cos((1-n)\theta) + i\sin((1-n)\theta))]

合并实部:

x=Re(zk+1)=n−1nrcos⁡θ+1nr1−ncos⁡((1−n)θ)x = \text{Re}(z_{k+1}) = \frac{n-1}{n}r\cos\theta + \frac{1}{n}r^{1-n}\cos((1-n)\theta)

⇒x=1n[(n−1)rcos⁡θ+r1−ncos⁡((1−n)θ)]\Rightarrow x=\frac{1}{n}[(n-1)r\cos \theta +r^{1-n}\cos\mathrm{((}1-n)\theta )]

合并虚部:

y=Im(zk+1)=n−1nrsin⁡θ+1nr1−nsin⁡((1−n)θ)y = \text{Im}(z_{k+1}) = \frac{n-1}{n}r\sin\theta + \frac{1}{n}r^{1-n}\sin((1-n)\theta)

⇒y=1n[(n−1)rsin⁡θ+r1−nsin⁡((1−n)θ)]\Rightarrow y=\frac{1}{n}[(n-1)r\sin \theta +r^{1-n}\sin\mathrm{((}1-n)\theta )]

即

zk=1n[(n−1)rcos⁡θ+r1−ncos⁡((1−n)θ)]+1n[(n−1)rsin⁡θ+r1−nsin⁡((1−n)θ)]iz_k=\frac{1}{n}[(n-1)r\cos \theta +r^{1-n}\cos\mathrm{((}1-n)\theta )]+\frac{1}{n}[(n-1)r\sin \theta +r^{1-n}\sin\mathrm{((}1-n)\theta )]i

‍

许可协议

本文采用 署名-非商业性使用-相同方式共享 4.0 国际 许可协议，转载请注明出处。

快来参与讨论吧~