隐式蒙特卡洛（Implicit Monte Carlo, IMC）方法最初由Fleck和Canfield在1971年提出，主要用于解决辐射传输问题。IMC方法通过结合传统的蒙特卡洛模拟与隐式时间积分的优点，允许更有效地模拟复杂的辐射输运过程。以下是一些科研领域中的具体应用。

1. 高能量密度物理

背景

高能量密度物理研究物质在极端条件下的行为，例如高温高压环境，这些条件常见于核爆炸、激光聚变实验等。

应用

在这种环境下，材料会经历快速的能量吸收和释放过程，导致温度和状态的急剧变化。IMC方法能够准确模拟辐射如何在这类环境中传播，并与物质发生相互作用，从而为设计实验提供理论支持。

2. 惯性约束聚变（ICF）

背景

惯性约束聚变是一种利用强激光或粒子束加热并压缩一个小燃料球（通常包含氘和氚），以达到核聚变条件的技术。

应用

在这个过程中，精确控制辐射场对于实现高效聚变至关重要。IMC方法可用于模拟不同阶段的辐射传输情况，帮助科学家优化激光脉冲形状、靶设计等关键参数，提高聚变效率。

3. 天体物理学

背景

天体物理学中常常需要研究恒星内部或星际介质中的辐射过程。

应用

IMC方法可以用来模拟恒星内部的能量传输机制，了解恒星演化的过程；也可以用于研究星际尘埃对背景辐射的吸收和再发射效应，增进我们对宇宙微波背景辐射的理解。

4. 核工程

背景

在核电站的设计和运行中，确保反应堆安全稳定地工作是至关重要的。

应用

IMC方法可以用于模拟反应堆内的辐射分布，评估不同材料对辐射的屏蔽效果，以及预测长期运行条件下结构材料的老化情况，从而指导新材料的选择和现有设施的安全评估。

参考文献

[1] Fleck Jr, J. A., & Canfield, T. R. (1971). Implicit Monte Carlo radiation transport in shock hydrodynamics. Journal of Computational Physics, 8(3), 313-342.

[2] Rosen, M. D. (1996). Inertial confinement fusion: Recent experimental and theoretical developments. Physics Today, 49(5), 40-45.

[3] Mihalas, D., & Weibel-Mihalas, B. (1984). Foundations of Radiation Hydrodynamics. Dover Publications.

[4] Duderstadt, J. J., & Hamilton, L. J. (1976). Nuclear Reactor Analysis. Wiley.

隐式蒙特卡洛方法通俗解释

隐式蒙特卡洛（Implicit Monte Carlo, IMC）方法是一种处理复杂系统中随机过程的高级数值模拟技术。为了帮助大家更好地理解IMC的核心思想，本文通过三个贴近生活的例子进行解释。

例子一：烤面包

想象你在做一个实验，目的是了解如何最有效地用烤箱烤好一片面包。在这个过程中，你有两个关键因素需要考虑：烤箱的温度和面包吸收热量并变热的速度。

传统方法的问题

如果按照传统的方法来做这个实验，你可能会频繁地检查面包的状态，并根据当前面包的热度调整烤箱的温度。这样做可能导致难以控制面包的加热过程，尤其是在快速加热而不让它烧焦的时候。

隐式蒙特卡洛方法

在每个步骤中，除了观察当前面包的热度外，你还尝试估计接下来一段时间内面包将会变得多么热。基于这种预估，你调整烤箱的温度设置或加热时间。比如，如果预估到面包很快就会太热，你就降低烤箱的温度；反之，如果觉得加热不够快，你就增加温度或延长加热时间。

例子二：加热游泳池

假设你有一个室外的大型游泳池，你想通过太阳能加热系统来加热水温。

传统方法的挑战

如果你使用传统的加热方法，你会根据当前水温来决定是否需要更多的热量。这种方法难以适应快速变化的环境条件，如天气突然变冷或变暖。

隐式蒙特卡洛方法的应用

首先，你不会只看当前的水温来做决策。相反，你会考虑未来的天气预报、一天中不同时间段的太阳强度等因素，预测接下来几个小时甚至几天内水温的变化趋势。基于这些预测，你可以更智能地调整太阳能加热系统的运行，以实现更加稳定和高效的加热过程。

例子三：和面

假设你需要准备2斤和好的面团，但是仅凭经验和直觉来加水和面粉（即传统的方法），可能会遇到以下问题：

传统方法

当你发现面团太湿时，就加面粉；如果面团太干，则加水。这样反复调整可能导致最终使用的材料远超过最初计划的量。

隐式蒙特卡洛方法

在开始和面前，你可以先做一个小规模的实验或根据经验预估，比如，如果你加了一定量的水后，估计需要多少面粉才能达到理想的湿度。基于上述预估，你可以做出更精确的调整，少量多次地添加面粉，并随时观察面团的状态，避免过度纠正带来的问题。

这三个例子共同展示了隐式蒙特卡洛方法的核心思想：不仅仅是基于当前状态做决策，还要考虑到未来可能的状态及其对当前决策的影响。这种方法允许我们在面对复杂系统时做出更加稳定和有效的调整，无论是物理模拟还是日常生活中的应用。

隐式（Implicit）这一概念在不同的领域有着不同的应用和含义，但总体来说，它指的是那些依赖于未来状态或间接定义的过程，而不是直接明确地给出解决方案。下面我们将从数值分析、编程及隐式蒙特卡洛方法三个角度来探讨隐式的含义。

数值分析中的隐式方法

定义与对比

在数值分析中，特别是在解决微分方程时，隐式方法是指那些在计算下一个时间步的解时需要用到该时间步未知解的方法。这意味着为了求得下一时刻的状态，我们需要通过求解一个方程组来获得结果，因为方程中包含了待求解的变量。

• 显式方法：基于当前已知的信息预测下一时间步的状态。例如，在向前欧拉法中，我们利用当前时刻的状态直接计算出下一时刻的状态。

• 隐式方法：考虑了未来时间点的信息。例如，在向后欧拉法中，为了确定下一时刻的状态，不仅需要当前时刻的状态信息，还需要该状态自身的未来值。这就导致了需要求解一个包含未知数的方程来找到确切的值。隐式方法的一个主要优点是可以使用较大的时间步长而不牺牲稳定性，特别适用于处理刚性问题。

编程中的隐式

在编程领域，“隐式”通常指不需要开发者直接指定所有细节的行为。这些行为由编译器、解释器或其他机制自动推断或执行。

• 隐式类型转换：比如在某些语言中，当两个不同类型的数值进行运算时，可能会自动将一种类型转换为另一种类型以完成操作，而无需程序员手动指定转换过程。例如，在JavaScript中，数字和字符串可以直接相加，其中数字会自动转换为字符串再进行连接。

• 隐式函数声明：在一些编程语言中，如果函数没有预先声明，编译器或解释器会根据函数调用的方式自动为其生成一个默认的声明。

隐式蒙特卡洛方法

隐式蒙特卡洛（Implicit Monte Carlo, IMC）是一种专门设计用来解决辐射传输问题的高级蒙特卡洛技术。IMC方法通过引入隐式时间离散化，显著提升了数值模拟的稳定性和准确性，尤其适用于处理复杂的辐射与物质相互作用问题。

• 传统蒙特卡洛方法：粒子（如光子）的轨迹是根据当前物理状态逐次随机抽样得到的。这种方法对于线性问题非常有效，但在处理非线性问题时可能遇到困难。

• 隐式蒙特卡洛方法：通过对粒子行为的概率分布进行调整，使得计算过程能够更准确地反映真实的物理过程。具体来说，IMC通过考虑未来时间点的状态信息，即在计算下一时刻的状态时，不仅依赖于当前已知的信息，还考虑到该状态自身的未来值，从而提高了计算的稳定性和准确性。

基本原理

蒙特卡洛方法求解定积分的理论依据是积分中值定理和大数定理。即

\[\int_a^bf(x)dx=f(\xi)\cdot(b-a)\]

其中\(a\leq\xi\leq b\). 蒙特卡洛方法是随机在\((a,b)\)区间取点，然后将这个点当作\(\xi\)来计算面积，但是由于点\(\xi\)当前是随机的，这个积分值就不是准确的，但是当大量取点时，由于函数\(f(x)\)本身的作用，根据大数定量最终会稳定到准确值的附近，因为接近真实值的情况出现的概率更大。

一维函数的定积分

假设我们想要计算如下函数在区间 \([0, 1]\) 上的积分:

\[\int_0^1x^2dx\]

该积分的解析解是\(\frac{1}{3}\). 使用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

import numpy as np

def monte_carlo_integral(f, a, b, N):
    """
    使用蒙特卡洛方法估计函数f在[a,b]上的积分。
    
    参数:
        f: 要积分的函数。
        a (float): 积分下限。
        b (float): 积分上限。
        N (int): 随机样本的数量。
    返回:
        float: 积分的蒙特卡洛估计值。
    """
    samples = np.random.uniform(a, b, N)  # 在[a, b]上生成N个均匀分布的随机数
    evaluations = f(samples)  # 计算每个随机样本点处函数的值
    integral_estimate = (b - a) * np.mean(evaluations)  # 积分估计
    return integral_estimate

# 定义要积分的函数
def f(x):
    return x**2

# 运行蒙特卡洛积分
a, b = 0, 1  # 积分范围
N = 100000   # 样本数量
integral = monte_carlo_integral(f, a, b, N)

print(f"蒙特卡洛积分估计值: {integral}")
print(f"真实值: {1/3}")

二维函数的定积分

二维函数的定积分同样也满足积分中值定理，所以蒙特卡洛方法是求解定积分的有效工具。本节我们求解定义在\([0,1]\times[0,1]\)上的二重积分:

\[\iint_{[0,1]\times[0,1]]}(x^2+y^2)dxdy\]

该积分的解可以直接解析求得为\(\frac{2}{3}\). 其 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

import numpy as np

def monte_carlo_double_integral(f, a1, b1, a2, b2, N):
    """
    使用蒙特卡洛方法估计函数f在[a1,b1] x [a2,b2]上的二重积分。
    
    参数:
        f: 要积分的函数，接受两个参数x和y。
        a1, b1 (float): 第一维的积分下限和上限。
        a2, b2 (float): 第二维的积分下限和上限。
        N (int): 随机样本的数量。
    返回:
        float: 积分的蒙特卡洛估计值。
    """
    # 在[a1, b1]和[a2, b2]上生成N个均匀分布的随机数对
    x_samples = np.random.uniform(a1, b1, N)
    y_samples = np.random.uniform(a2, b2, N)
    
    # 计算每个随机样本点处函数的值
    evaluations = f(x_samples, y_samples)
    
    # 积分估计
    area = (b1 - a1) * (b2 - a2)
    integral_estimate = area * np.mean(evaluations)
    
    return integral_estimate

# 定义要积分的函数
def integrand(x, y):
    return x**2 + y**2

# 运行蒙特卡洛二重积分
a1, b1 = 0, 1  # 第一维积分范围
a2, b2 = 0, 1  # 第二维积分范围
N = 100000     # 样本数量

integral = monte_carlo_double_integral(integrand, a1, b1, a2, b2, N)

print(f"蒙特卡洛积分估计值: {integral}")
print(f"真实值: {2/3}")

使用蒙特卡洛方法求解一阶线性常微分方程

在这个例子中，我们将展示如何使用蒙特卡洛方法来近似求解一个简单的一阶线性常微分方程，并与解析解进行比较。

问题描述

考虑如下形式的初值问题：

\[\frac{dy}{dx}=-2xy,\quad y(0)= 1\]

其解析解为 \(y(x) = e^{-x^2}\)。我们将使用蒙特卡洛方法来近似这个解。

Python代码实现

以下是基于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

import numpy as np
import matplotlib.pyplot as plt

def monte_carlo_ode_solver(f, x0, y0, x_end, num_samples=1000, dx=0.05):
    """
    使用蒙特卡洛方法求解一阶常微分方程。
    
    参数:
        f: 函数，定义了dy/dx = f(x, y)
        x0 (float): 初始x值。
        y0 (float): 初始y值。
        x_end (float): x的最大值。
        num_samples (int): 每一步中使用的样本数量。
        dx (float): x方向上的步长。
    返回:
        x_points (numpy.ndarray): x坐标点。
        y_estimated (numpy.ndarray): 对应每个x点估计的y值。
    """
    x_points = np.arange(x0, x_end+dx, dx)
    y_estimated = np.zeros_like(x_points)
    y_estimated[0] = y0  # 设置初始条件
    
    for i in range(1, len(x_points)):
        current_x = x_points[i-1]
        current_y = y_estimated[i-1]
        
        dy_samples = []
        for _ in range(num_samples):
            # 随机扰动dx
            delta_x = dx * (np.random.rand() - 0.5)
            # 根据当前点的斜率计算dy
            delta_y = f(current_x, current_y) * (dx + delta_x)
            dy_samples.append(delta_y)
        
        # 计算平均dy并更新y值
        avg_dy = np.mean(dy_samples)
        y_estimated[i] = current_y + avg_dy
    
    return x_points, y_estimated

# 定义导数函数
def f(x, y):
    return -2 * x * y

# 运行求解器
x, y_mc = monte_carlo_ode_solver(f, x0=0, y0=1, x_end=2, num_samples=1000)

# 计算真实解用于比较
y_true = np.exp(-x**2)

# 绘制结果
plt.plot(x, y_mc, label='Monte Carlo Approximation')
plt.plot(x, y_true, label='True Solution', linestyle='--')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show() 

分析

在使用蒙特卡洛方法求解微分方程时，最关键的一步是 delta_x = dx * (np.random.rand() - 0.5), 这表示对于步长\(dx\)在其中心附近扰动，并且保存每次扰动的所造成的\(dy\)值，当完成大量的数据采样时，根据大数定理，\(dy\)值会稳定到真实的\(dy\)的情况。进一步将此值追加到\(y\)，同时\(dx\)追加到\(x\)后，再进行下一个点的判断。所以蒙特卡洛方法基于大数定理是一个有效的方法。

使用蒙特卡洛方法估算圆周率

在这篇博客文章中，我们将探讨如何使用蒙特卡洛方法来估算圆周率 \(\pi\) 的值，并提供一个简单的Python代码示例以及详细的代码解析。

蒙特卡洛方法简介

蒙特卡洛方法是一种通过随机抽样进行数值计算的方法。它广泛应用于物理、工程、金融等领域。本文将演示如何利用这种方法估算圆周率 \(\pi\)。

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

# 导入random模块，用于生成伪随机数
import random

# 定义函数estimate_pi，接受参数num_samples表示要生成的随机点的数量
def estimate_pi(num_samples):
    # 初始化计数器变量num_points_in_circle为0，用于记录落在单位圆内的点的数量
    num_points_in_circle = 0
    
    # 使用for循环迭代num_samples次，下划线_表示我们不关心循环变量的具体值
    for _ in range(num_samples):
        # 在每次循环中，使用random.uniform(-1, 1)生成两个介于-1到1之间的随机浮点数，分别赋给变量x和y
        x, y = random.uniform(-1, 1), random.uniform(-1, 1)
        # 如果该点位于半径为1的单位圆内（即检查x^2 + y^2 <= 1），则增加num_points_in_circle计数器的值
        if x**2 + y**2 <= 1:
            num_points_in_circle += 1
    
    # 根据公式4 * num_points_in_circle / num_samples来估算π的值
    pi_estimate = 4 * num_points_in_circle / num_samples
    # 返回计算出的π的估计值
    return pi_estimate

# 设置模拟次数为100,000
num_samples = 100000
# 调用estimate_pi函数获取π的估计值
pi_estimate = estimate_pi(num_samples)

# 输出结果，显示使用了多少个样本以及基于这些样本估计得到的π的值
print(f"Estimated value of Pi using {num_samples} samples: {pi_estimate}") 

当用户从Windows 10/11切换到Arch Linux时，如果发现无线网卡不见了，这通常是由于Windows的快速启动特性导致的。快速启动会使得某些硬件（如无线网卡）处于一种休眠状态，这样当你切换到Linux系统时，这些硬件可能不会被正确唤醒。

解决方案

• 对于双系统的用户：建议去Windows那边执行一次真正的重启（而不是简单地关机再开机），因为Windows的“关机”选项实际上可能是快速启动，它并不会完全关闭计算机而是让系统进入一个类似休眠的状态。通过在Windows中选择“重启”，可以确保所有设备都被正确重置，包括无线网卡。

  • 进入Windows，然后点击“重启”而非“关机”。这通常能解决无线网卡在Linux中不显示的问题。

• 对于只有Arch Linux的用户：如果您的机器只安装了Arch Linux并且也遇到了同样的问题，那么问题的原因可能与其他因素有关，比如驱动程序、内核版本或者是电源管理设置等。需要进一步检查系统日志或尝试更新系统及固件。

• 一劳永逸的解决方案：虽然通过重启Windows可以暂时解决问题，但为了从根本上避免这种情况的发生，您可以考虑在Windows中禁用快速启动功能。这样每次从Windows切换到Linux时都不会受到快速启动的影响，从而避免无线网卡无法被识别的问题。

  • 在管理员命令提示符窗口，输入：powercfg /h on命令,启用Windows快速启动；
  • 输入：powercfg /h off命令, 关闭Windows快速启动。

参考文章

• Arch Linux中文论坛讨论
• 如何关闭或禁用Windows 10或11的快速启动！
• Windows10系统启用和关闭快速启动的两种方法

Mandelbrot 集简介

Mandelbrot 集是数学中一个非常著名的分形集合，它以法国数学家本华·曼德勃罗的名字命名。分形是一种具有自相似性质的几何形状，意味着无论放大多少倍，其部分看起来与整体相似。 Mandelbrot 集正是这种特性的典型代表，并且因其复杂而美丽的图形成为了计算机图形学、艺术以及数学研究中的重要主题。

Mandelbrot 集定义于复数平面上，由所有不使序列 $ z_{n+1} = z_n^2 + c $ 发散至无穷大的点 $ c $ 组成，其中 $ z_0 = 0 $。换句话说，对于每一个位于复数平面上的点 $ c $，从 $ z=0 $ 开始应用上述迭代公式，如果该序列保持有界，则点 $ c $ 属于 Mandelbrot 集；反之，如果序列趋向于无穷大，则该点不属于此集合。通常情况下，为了判断一个点是否属于 Mandelbrot 集，会设定一个最大迭代次数和逃逸半径（通常是2），如果在达到最大迭代次数前，序列的模超过了逃逸半径，则认为该点发散。

直接实现方式

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

#! /usr/bin/env python3
# vim:fenc=utf-8
import numpy as np
# Mandelbrot Set参数
picname = "Mandelbrot"
file_name = f"{picname}.csv"
width, height = 1000, 1000  # 图像的宽度和高度
max_iter = 100  # 最大迭代次数

def mandelbrot_set(width, height, max_iter):
    x = np.linspace(-2.5, 1.5, width)
    y = np.linspace(-2, 2, height)
    escape_time = np.zeros((height, width))  # 注意这里的维度顺序

    for i in range(width):
        for j in range(height):
            c = complex(x[i], y[j])
            z = 0
            n = 0
            while abs(z) <= 2 and n < max_iter:
                z = z*z + c
                n += 1
            escape_time[j, i] = n  # 根据像素位置更新逃逸时间
            
    return x, y, escape_time

# 计算Mandelbrot集
x_coords, y_coords, escape_time = mandelbrot_set(width, height, max_iter)

# 将结果保存为适合Veusz导入的格式
with open(file_name, 'w') as file:
    # 写入X坐标行
    file.write(',' + ','.join(map(str, x_coords)) + '\n')
    for j in range(escape_time.shape[0]):
        # 写入每行的数据，先写Y坐标，然后是该行的逃逸时间值
        file.write(str(y_coords[j]) + ',' + ','.join(map(str, escape_time[j])) + '\n')

print(f"Mandelbrot集已成功保存到{file_name}")

使用 Pandas 库保存数据

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

#! /usr/bin/env python3
# vim:fenc=utf-8
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

def mandelbrot(c, max_iter):
    z = 0
    n = 0
    while abs(z) <= 2 and n < max_iter:
        z = z*z + c
        n += 1
    return n

def mandelbrot_set(xmin, xmax, ymin, ymax, width, height, max_iter):
    r1 = np.linspace(xmin, xmax, width)
    r2 = np.linspace(ymin, ymax, height)
    m = np.empty((width,height))
    for i in range(width):
        for j in range(height):
            m[i,j] = mandelbrot(r1[i] + 1j*r2[j], max_iter)
    return m, r1, r2  # 返回m以及对应的x和y坐标值

# 设置参数并生成Mandelbrot集合
xmin, xmax, ymin, ymax = -2.0, 1.0, -1.5, 1.5
width, height = 800, 800
max_iter = 256

m, r1, r2 = mandelbrot_set(xmin, xmax, ymin, ymax, width, height, max_iter)

# 显示图像（可选）
# plt.imshow(m.T, extent=[xmin, xmax, ymin, ymax], cmap='hot', interpolation='bilinear')
# plt.colorbar()
# plt.title("Mandelbrot Set")
# plt.show()

# 将结果保存为CSV文件
# 创建一个列表来保存数据点
data_points = []
for i in range(width):
    for j in range(height):
        data_points.append([r1[i], r2[j], m[i, j]])

# 使用pandas DataFrame保存数据
df = pd.DataFrame(data_points, columns=['Re', 'Im', 'Iteration'])
df.to_csv('mandelbrot_set.csv', index=False)

print("Mandelbrot集合已成功保存至mandelbrot_set.csv")

Pandas 保存数据与直接保存数据的不同

在处理和存储数据时，使用Pandas的DataFrame来保存数据与直接保存数据（例如将数据保存为纯文本、CSV文件等）有着显著的区别。以下是这两种方法的主要不同点：

使用 Pandas DataFrame保存数据

结构化存储

• 表格形式：DataFrame提供了一个二维表格结构，支持行和列标签，并允许每列有不同的数据类型。
• 强大的数据操作能力：内置了丰富的函数用于数据清洗、转换、筛选、聚合等高级功能。

索引和标签支持

• 灵活访问：可以通过列名、行位置或布尔索引来选取子集，便于快速定位所需信息。
• 多级索引：支持复杂的索引机制，使得数据管理和查询更加灵活。

缺失值处理

• 自动管理缺失数据：通过NaN表示缺失值，并提供了多种方法来清理或填补这些缺失值。

数据可视化和分析集成

• 与其他库的良好集成：可以方便地与NumPy、Matplotlib、Seaborn等科学计算和可视化库结合使用，增强数据分析能力。

直接保存数据

简单的数据格式

• 通用性高：直接存储的数据通常保存为原始格式，如CSV、TXT、JSON等，适合于结构相对简单、规模较小的数据集。
• 轻量级：这种保存方式不依赖于外部库，适合需要跨平台使用的场景。

有限的数据处理能力

• 缺乏高级特性：直接存储的数据没有内置的数据处理功能，如果需要进行数据处理，则需编写额外的代码或使用其他工具。
• 基础操作：虽然可以直接读写数据，但复杂的数据操作如过滤、分组、聚合等需要额外的努力来实现。

通用性和可移植性

• 易于分享：由于直接存储的数据格式通常是通用的，因此可以在不同的编程语言和环境中读取和写入，增加了跨平台使用的便利性。

总的来说，选择哪种方式取决于具体的应用需求和数据特点。DataFrame 提供了一种更为强大和灵活的方式来存储和处理数据，特别适用于复杂的数据分析任务；而直接保存数据更适合于简单场景，或者当你需要一种轻量级、跨平台的数据交换方式时。

Julia 集简介

朱利亚集是分形几何中的一个重要概念，由法国数学家加斯东·朱利亚研究得出，通过迭代复数平面上的函数 \(f_c(z) = z^2 + c\) 来定义，其中 \(c\) 为固定的复数参数。对于每一个给定的 \(c\) 值，从任意初始点 \(z_0\) 开始进行迭代，如果序列 \({zn}\) 保持有界，则初始点属于该朱利亚集；反之，若序列趋向无穷大，则初始点不属于此集。这些集合并呈现出高度复杂和自相似的边界，其形态随不同的 \(c\) 值变化，从连通区域到复杂的尘埃状分布不等，展现了丰富的几何结构和美丽的图案，既是数学研究的对象，也是艺术创作的源泉。

通过计算机绘制 Julia 集，需要指定一个 \(c\) 值， 然后取一个复平面上的一个区域 \(A\)， 对于此区域内的每一点按\(f(z)\)进行迭代, 同时设置一个最大迭代次数 \(N\) 和逃逸半径 \(R\). 对于复平面内的点如果\(|z|>R\), 则取强度为\(0\); 如果 \(|z|<R\), 则使用 \(f(z)\) 迭代， 当迭代 \(m\) 次后点 \(|z_m|>R\), 则取强度为\(m\), 如果迭代\(N\)次后仍然不满足\(|z_N|>R\), 则取强度为\(N\). 如此操作后，对于复平面区域\(A\)内的每一点都赋予一个强度值，然后输出为 csv 的2D数据，按强度上色就构成了一个 Julia 集的图形显示。

代码分析

声明脚本为 python3 , 用 utf-8 编码

1
2

#! /usr/bin/env python3
# vim:fenc=utf-8

调用 numpy 库, 同时设置初值 c, 并定义保存文件的名称 picname

1
2
3
4
5
6

import numpy as np
c = complex(-0.70176, -0.3842)
picname = "Nebula"
file_name = f"{picname}.csv"
width, height = 1000, 1000  # 图像的宽度和高度
max_iter = 100  # 最大迭代次数

绘制程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14

def julia_set(c, width, height, max_iter):
    x = np.linspace(-2, 2, width)
    y = np.linspace(-2, 2, height)
    escape_time = np.zeros((height, width))  # 注意这里的维度顺序
    for i in range(width):
        for j in range(height):
            z = complex(x[i], y[j])
            n = 0
            while abs(z) <= 2 and n < max_iter:
                z = z*z + c
                n += 1
            escape_time[j, i] = n  # 根据像素位置更新逃逸时间
            
    return x, y, escape_time

上述程序中：

• 第1行：定义程序julia_set, 并设置了四个参数c, width, height 和 max_iter.
• 第2行：定义x数组，将x定义域\(-2\)到\(2\)划分成width等份.
• 第3行：定义y数组，将y定义域\(-2\)到\(2\)划分成height等份.
• 第4行：定义 \(height\times width\) 矩阵escape_time，所有矩阵元默认取 \(0\).
• 第5行：遍历x数组, 以取得每一个 \(x\)数组元 \(x[i]\).
• 第6行：遍历y数组, 以取得每一个 \(y\)数组元 \(y[j]\).
• 第7行：以\(x[i]\) 为实部， \(y[j]\)为虚部，构建复数\(z\).
• 第8行：对于每一个点\(z\), 默认逃逸次数为\(n=0\).
• 第9行：使用while对每一个点\(z\)判断，条件为小于逃逸半径2和小于最大逃逸次数 max_iter.
• 第10行：若满足条件，则迭代一次，同时逃逸次数 \(n\) 增加\(1\), 若仍然满足第9行的条件，则继续循环，同时不断将判断的逃逸次数\(n\)赋值给矩阵escape_time[j,i].
• 第11行：返回最后的x,y,escape_time 数组.

使用上述函数处理具体的参数

1
2

# 计算Julia集
x_coords, y_coords, escape_time = julia_set(c, width, height, max_iter)

将结果保存为 csv 数据

1
2
3
4
5
6

with open(file_name, 'w') as file:
    file.write(',' + ','.join(map(str, x_coords)) + '\n')
    for j in range(escape_time.shape[0]):
        # 写入每行的数据，先写Y坐标，然后是该行的逃逸时间值
        file.write(str(y_coords[j]) + ',' + ','.join(map(str, escape_time[j])) + '\n')
print(f"Julia集已成功保存到{file_name}")

上述程序中：

• 第1行：打开文件file_name, 同时赋予write权限（即选项w）, 将其命名为file， 在后续程序中调用.
• 第2行：写入x坐标行, 数组以逗号,分隔, 行尾添加换行符\n.
• 第3行：遍历escape_time数组的第一列（也就是所有行数）, 写入每一行的数组，仍然以逗号,分隔，行尾添加换行符\n.
• 第4行：打印处理结果，同时使用f标记来使用变量{file_name}.

今天需要使用 veusz 处理一个 julia 集的二维图形，但是发现 veusz 不能正常启动了。由于 veusz 是一个使用 python 写的软件， 所以由于 archlinux 升级导致的问题可能性还是比较大的, 于是卸载 veusz 后重新安装，得到了错误提示：

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

Traceback (most recent call last):
  File "/usr/lib/python3.13/site-packages/pyproject_hooks/_in_process/_in_process.py", line 389, in <module>
    main()
    ~~~~^^
  File "/usr/lib/python3.13/site-packages/pyproject_hooks/_in_process/_in_process.py", line 373, in main
    json_out["return_val"] = hook(**hook_input["kwargs"])
                             ~~~~^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3.13/site-packages/pyproject_hooks/_in_process/_in_process.py", line 143, in get_requires_for_build_wheel
    return hook(config_settings)
  File "/usr/lib/python3.13/site-packages/setuptools/build_meta.py", line 332, in get_requires_for_build_wheel
    return self._get_build_requires(config_settings, requirements=[])
           ~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3.13/site-packages/setuptools/build_meta.py", line 302, in _get_build_requires
    self.run_setup()
    ~~~~~~~~~~~~~~^^
  File "/usr/lib/python3.13/site-packages/setuptools/build_meta.py", line 516, in run_setup
    super().run_setup(setup_script=setup_script)
    ~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3.13/site-packages/setuptools/build_meta.py", line 318, in run_setup
    exec(code, locals())
    ~~~~^^^^^^^^^^^^^^^^
  File "<string>", line 38, in <module>
  File "/home/feng/.cache/paru/clone/veusz/src/veusz-3.6.2/pyqt_setuptools.py", line 11, in <module>
    import tomli
ModuleNotFoundError: No module named 'tomli'

显然提示缺少 python 库 tomli, 于是使用 archlinux 自带的包管理器 pacman 安装 tomli

1

sudo pacman -S python-tomli

之后再安装 veusz, 则问题解决。