BNU-FZH

fengzhenhua@outlook.com

在使用蒙特卡洛方法计算圆周率时,选择积分公式的关键在于方差的大小实现的效率。本文通过数学推导和实验验证,对比两种常见积分公式的优劣。


1. 积分公式的数学等价性

两个积分的数学结果均为π,但被积函数特性不同:

  • 公式1
    \[ 2 \int_0^1 \frac{1}{\sqrt{1 - x^2}} \, dx = \pi \]
    被积函数 \(f_1(x) = \frac{1}{\sqrt{1 - x^2}}\)\(x \to 1\) 时发散,导致数值计算的不稳定性。

  • 公式2
    \[ 2 \int_{-1}^1 \sqrt{1 - x^2} \, dx = \pi \]
    被积函数 \(f_2(x) = \sqrt{1 - x^2}\) 在区间 \([-1, 1]\) 上有界且连续,方差更小。


2. 蒙特卡洛方法的方差分析

2.1 公式1的方差问题

  • 发散特性
    \(f_1(x)\)\(x \to 1\) 处趋向无穷大,导致蒙特卡洛采样时方差显著增大。例如,当 \(x = 0.99\) 时,\(f_1(x) \approx 7\),远高于平均值 \(\pi/2 \approx 1.57\),极端值会严重干扰估计结果。

  • 实现挑战
    需要均匀采样 \(x \in [0, 1]\),但靠近 \(x = 1\) 的点对结果贡献极大,易引入噪声。

2.2 公式2的方差优势

  • 有界性
    \(f_2(x)\) 的值域为 \([0, 1]\),方差远小于公式1。例如,\(N=10^6\) 次采样时,公式2的估计误差通常比公式1低一个数量级。

  • 数值稳定性
    无需处理函数发散问题,所有采样点的贡献均在可控范围内。


3. 实验验证:方差的量化对比

假设使用 \(N\) 次采样:

  • 公式1
    由于 \(f_1(x)\) 的发散性,方差 \(D_1 \propto \int_0^1 \left( \frac{1}{\sqrt{1-x^2}} \right)^2 dx = \infty\)(发散)。

  • 公式2
    方差 \(D_2 = \int_{-1}^1 \left( \sqrt{1-x^2} \right)^2 dx - \left( \frac{\pi}{2} \right)^2 = \frac{\pi}{2} - \frac{\pi^2}{4} \approx 0.43\),远小于公式1。


4. 代码实现对比

4.1 公式2的C++实现(高效且低方差)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import random
import math

def estimate_pi(n):
count = 0
for _ in range(n):
x = random.uniform(-1, 1) # 生成x ∈ [-1, 1]
y = random.uniform(0, 1) # 生成y ∈ [0, 1]
if y <= math.sqrt(1 - x**2): # 判断是否在半圆内
count += 1
# 计算公式:2 * (积分值) = 2 * (count/N * 矩形面积)
return 2.0 * (count * 2.0 / n) # 矩形面积为2(宽2×高1)

# 示例:使用100万个样本
print(f"Estimated π: {estimate_pi(1000000)}")

4.2 公式1的Python实现(高方差问题)

1
2
3
4
5
6
7
8
9
10
import numpy as np

def monte_carlo_formula1(n):
x = np.random.uniform(0, 1, n)
f = 1 / np.sqrt(1 - x**2)
integral = np.mean(f)
return 2 * integral

# 测试示例(方差极大)
print(monte_carlo_formula1(1000000)) # 结果可能波动明显

5. 结论

公式2(\(2\int_{-1}^{1} \sqrt{1 - x^2} \, dx\))是更优选择,原因如下:

  1. 方差更小:被积函数 \(\sqrt{1 - x^2}\) 在区间 \([-1, 1]\) 上有界且连续,使得蒙特卡洛方法的收敛速度更快。
  2. 数值稳定性:无需处理函数发散问题(如公式1中 \(\frac{1}{\sqrt{1-x^2}}\)\(x \to 1\) 处的发散),避免极端值对估计结果的干扰。
  3. 计算效率:在相同采样次数下,公式2的估计误差显著低于公式1,且计算过程更稳定。

选择积分公式时,应优先考虑被积函数的有界性光滑性,以降低蒙特卡洛方法的方差,从而提升计算效率。

如果 TeXstudio 出现问题无法打开或闪退,你可以尝试恢复其默认设置。TeXstudio 将其设置存储在一个名为 texstudio.ini 的配置文件中。要恢复默认设置,你需要删除或重命名这个文件,这样 TeXstudio 在下次启动时会生成一个新的默认配置文件。

操作步骤

1. 关闭 TeXstudio

确保 TeXstudio 已经完全关闭。如果它因为错误而无法正常关闭,请通过任务管理器强制结束它的进程。

2. 找到配置文件

  • Windows: 通常位于用户目录下的 AppData\Roaming\texstudio\texstudio.ini。注意 AppData 文件夹是隐藏的,你可能需要在资源管理器中显示隐藏的文件和文件夹才能看到它。
  • macOS: 配置文件一般位于 ~/Library/Preferences/texstudio.ini
  • Linux: 对于大多数发行版来说,配置文件可以在 ~/.config/texstudio/texstudio.ini 找到。

3. 删除或重命名配置文件

为了安全起见,你可以先尝试重命名该文件(例如,添加 .bak 后缀),而不是直接删除它。这样,如果你发现新的默认设置有问题,还可以恢复原来的配置。

4. 重新启动 TeXstudio

现在再次启动 TeXstudio,程序将会使用默认设置启动。

5. 验证更改

检查是否所有的问题都得到了解决,并根据需要重新配置你的偏好设置。

注意:此过程将删除所有自定义设置,包括快捷键、语法高亮设置等。如果你有重要的自定义设置,建议在执行上述步骤前备份相关的配置信息。

LibreOffice是一款开源的办公软件套件,能够提供文字处理、电子表格、演示文稿、绘图、数据库等多种功能,适用于Windows、MacOS以及Linux等多种操作系统。然而我在ArchLinux使用LibreOffice时遇到一些问题,本文记录这些问题的解决方案。

取消表格选中时高对比度

2025年03月14日星期五阴北京市, 当安装好LibreOffice, 默认开启了高对比度,这导致每次选中表格的单元格时都是黑色的,相当不雅观。解决方法为:工具选项LibreOffice无障碍辅助高对比度禁用.

工具栏图标模糊的解决方法

2025年03月14日星期五阴北京市, 在4k显示器使用ArchLinux系统,LibreOffice默认开启的工具栏图标不是svg格式,这导致了工具栏图标模糊,解决方法为:工具选项LibreOffice视图主题选择任何一款svg图标.

在泰勒展开中,余项的形式决定了展开式的精度和适用范围。拉格朗日余项(Lagrange Remainder)和皮亚诺余项(Peano Remainder)是两种最常见的余项类型,它们的核心区别在于对误差的描述方式和无穷小性的要求。本文将通过表达式、性质、对比和例子详细分析两者的异同。


1. 拉格朗日余项(Lagrange Remainder)

表达式

对于函数 $ y(x) $ 在 $ x_0 $ 处的 $ n $ 阶泰勒展开,拉格朗日余项的形式为: \[ R_n(x) = \frac{y^{(n+1)}(\xi)}{(n+1)!}(x - x_0)^{n+1} \] 其中,$ $ 是介于 $ x_0 $ 和 $ x $ 之间的某个点。

性质

  • 余项的无穷小性: 当 $ x_0$ 时,若 $ y^{(n+1)}(x) $ 在 $ x_0 $ 附近有界,则余项 $ R_n(x) $ 趋于零;否则余项可能不趋于零。
  • 误差描述: 提供定量的误差估计,例如明确给出余项的表达式和依赖的高阶导数值。
  • 导数要求: 需要函数 $ y(x) $ 在 $ x_0 $ 处具有 $ n+1 $ 阶连续导数

例子

以 $ y(x) = e^x $ 在 $ x_0 = 0 $ 处的二阶泰勒展开为例: \[ e^x = 1 + x + \frac{x^2}{2} + \frac{e^\xi}{6}x^3 \quad (\xi \text{ 介于 } 0 \text{ 和 } x \text{ 之间}) \] 余项 $ x^3 $ 的大小取决于 $ $ 的值。当 $ x $ 时,$ $,余项趋于零;但对固定 $ x $,余项不一定无穷小。


2. 皮亚诺余项(Peano Remainder)

表达式

对于 $ n $ 阶泰勒展开,皮亚诺余项的形式为: \[ R_n(x) = o\left((x - x_0)^n\right) \quad \text{当} \ x \to x_0 \] 符号 $ o(...) $ 表示余项是比 $ (x - x_0)^n $ 更高阶的无穷小

性质

  • 余项的无穷小性
    强制要求余项满足: \[ \lim_{x \to x_0} \frac{R_n(x)}{(x - x_0)^n} = 0 \] 即余项必须比展开式的最后一项 $ (x - x_0)^n $ 更快趋近于零。
  • 误差描述
    仅提供定性描述,强调余项的高阶无穷小性质,不涉及具体数值估计。
  • 导数要求
    仅需函数 $ y(x) $ 在 $ x_0 $ 处存在 ​$ n $ 阶导数,无需更高阶导数连续。

例子

以 $ y(x) = e^x $ 在 $ x_0 = 0 $ 处的二阶泰勒展开为例: \[ e^x = 1 + x + \frac{x^2}{2} + o(x^2) \quad (x \to 0) \] 余项 $ o(x^2) $ 仅表明当 \(x \to 0\) 时,余项比 $ x^2 $ 更快趋近于零,但不提供具体表达式。


关键区别对比

性质 拉格朗日余项 皮亚诺余项
余项是否无穷小 不一定,取决于高阶导数的有界性 必须无穷小(定义要求)
误差描述 定量估计(具体数值) 定性描述(仅高阶无穷小)
导数要求 需要 $ y^{(n+1)}(x) $ 存在且连续 只需 $ y^{(n)}(x) $ 在 $ x_0 $ 处存在
适用场景 需要误差范围的实际问题 理论分析,强调局部逼近性质

结论

  1. 拉格朗日余项
    • 适用于需要定量误差估计的场景(如工程计算)。
    • 余项的无穷小性依赖于高阶导数的有界性,当 \(x \to x_0\) 且导数有界时,余项趋于零。
  2. 皮亚诺余项
    • 适用于理论分析,强调函数在局部的高阶逼近性质。
    • 强制余项为高阶无穷小,无需具体计算高阶导数的值。

两种余项的共同点是当 \(x \to x_0\) 时,余项均趋于零,但皮亚诺余项对无穷小性的要求更严格,且不依赖高阶导数的连续性。


延伸思考: 在实际应用中,若函数的高阶导数难以计算或无界(例如 $ y(x) = x^{5/2} $ 在 $ x=0 $ 处的三阶导数不存在),皮亚诺余项可能是唯一可行的选择。

Gtk4 Desktop Icons NG (DING) 简介

Gtk4 Desktop Icons NG (DING) 是一个GNOME Shell扩展,旨在为使用GTK4的GNOME桌面环境提供增强的桌面图标管理功能。这个扩展是对原有支持GTK3的Desktop Icons NG的一个更新版本,它带来了对GTK4的支持,并引入了一系列改进和新特性,以提升用户体验。

Gtk4 Desktop Icons NG

主要特点包括:

  • 灵活的图标布局:图标可以自由放置在桌面上任何位置,或者选择对齐到网格,方便用户根据个人喜好组织桌面。

  • 桌面链接创建:允许用户直接在桌面上创建文件或文件夹的快捷方式,便于快速访问重要资源。

  • GSconnect集成:通过与GSconnect的集成,用户可以直接从桌面发送文件到连接的设备(如手机),极大地提高了跨设备文件传输的便捷性。

  • 拖放支持:支持将文件拖放到Dock、Dash上,或反向操作,即从Dock、Dash拖放到桌面上,增强了操作灵活性。

  • 现代化的代码库:采用了更新的技术栈,包括移植到ESM模块和使用Gio菜单,确保了更好的性能和更高的稳定性。所有翻译工作都在Weblate平台上进行,确保多语言支持的质量。

  • 异步函数:尽可能采用异步处理的方式实现各种功能,提升了响应速度和整体性能。

  • 兼容性:支持Gnome 45及更高版本,确保用户可以在最新的GNOME环境中享受这些改进。

总的来说,Gtk4 Desktop Icons NG (DING) 提供了一个更加现代化、功能丰富的解决方案来管理和定制你的GNOME桌面环境中的图标显示,特别适合那些希望获得更传统桌面体验并享受最新技术带来的优势的用户。

插件安装

  1. 访问插件网站:GNOME Extensions.
  2. 搜索Gtk4 Desktop Icons NG (DING), 点击右侧的按钮On.

使用蒙特卡洛方法求解圆周率的一个经典方法是取一个单位圆和外接正四边形,则正四边形的面积为1, 随机生成点数,统计落在圆内的点的数目m, 然后再以m除以总点数n,获得随机点落在圆内的概率,再乘以正方形的面积就得到了圆周率。 但是本文不计划采用这个方法,我们使用积分来计算,即

\[\pi=2\int_0^1 \frac{1}{\sqrt{1-x^2}} dx\]

按积分中值定理,在\((0,1)\)内必存在一点\(\xi\)满足

\[\pi=2\frac{1}{\sqrt{1-\xi^2}}(1-0)\]

如果我们任意取\((0,1)\)中的一个点,那必然不一定满足上式,但是按照大数定理,当在\((0,1)\)内随机取值的数量达到足够大,并且平均时,按统计规律积分值将接近真实结果,并且随随机取点的数量增大而提高精度。

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
#! /usr/bin/env python3
# vim:fenc=utf-8
import numpy as np
# 定义函数f(x)
def f(x):
return 1/(1-x**2)**0.5

# Monte Carlo方法求积分
def monte_carlo_integral(f, a, b, num_samples=90000000):
# 在区间[a, b]内随机生成num_samples个样本点
samples = np.random.uniform(a, b, num_samples)

# 计算这些样本点对应的函数值
y_values = f(samples)

# 根据平均高度估计积分值,乘以区间的宽度得到近似面积
integral_estimate = (b - a) * np.mean(y_values)

return 2*integral_estimate

# 使用Monte Carlo方法估计积分
a, b = 0, 1 # 定义积分区间
estimated_integral = monte_carlo_integral(f, a, b)

print("Pi", estimated_integral)

Rust 实现

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
use rand::distributions::{Distribution, Uniform};
use std::f64;

fn f(x: f64) -> f64 {
1.0 / (1.0 - x.powi(2)).sqrt()
}

fn monte_carlo_integral(a: f64, b: f64, num_samples: usize) -> f64 {
let between = Uniform::new(a, b);
let mut rng = rand::thread_rng();

let mut sum = 0.0;
for _ in 0..num_samples {
let sample = between.sample(&mut rng);
sum += f(sample);
}

let integral_estimate = (b - a) * (sum / num_samples as f64);
2.0 * integral_estimate
}

fn main() {
let a = 0.0;
let b = 1.0;
let num_samples = 90_000_000;

let estimated_integral = monte_carlo_integral(a, b, num_samples);

println!("Pi: {}", estimated_integral);
}

在这个例子中,我们使用了 rand crate 来生成随机数。因此,在你的 Cargo.toml 文件中需要添加以下依赖:

1
2
[dependencies]
rand = "0.8"

使用cargo build --release编译后再运行代码,可以明显感觉到比Python要快的多。

  • 运行环境:ArchLinux,ThinkPad T490, Intel(R) Core(TM) i5-8265U (8) @ 3.90 GHz,NVIDIA GeForce MX250 [Discrete],Memory 16G
  • Python 执行总用时1.255s
  • Rust 执行总用时0.466s
  • 它们的精度都能达到3.14xxx, 通过这个例子也显示出了Rust的优越性,在这个例子它需要的时间只是Python的0.37倍。

理解 Rust 的所有权系统

Rust 的所有权(Ownership)是其核心特性之一,旨在确保内存安全的同时提供高性能。通过所有权系统,Rust 能够在编译时自动管理内存,避免了诸如空指针引用、数据竞争等常见问题,而无需垃圾回收机制。

什么是所有权?

每个程序都需要管理它使用的内存资源。在低级语言中,如 C 或 C++,程序员必须手动分配和释放内存,这容易导致错误。Rust 使用所有权系统来自动化这一过程,并且保证了内存的安全性。

主要概念

  • 所有权:控制着堆上数据的生命周期。
  • 借用:允许临时访问某个值而不获取其所有权。
  • 生命周期:定义了引用的有效范围,确保引用始终有效。

所有权规则

  1. 每个值都有一个变量,称为它的所有者。
  2. 每个值在同一时间只能有一个所有者。
  3. 当所有者离开作用域时,该值将被丢弃。

示例代码

1
2
3
4
5
6
7
fn main() {
let s1 = String::from("hello"); // s1 是所有者
let s2 = s1; // s1 的所有权转移给 s2,s1 不再有效

// println!("{}", s1); // 这行会报错,因为 s1 已经不再有效
println!("{}", s2); // 正确,s2 拥有 "hello"
}

借用

为了避免所有权转移带来的不便,Rust 允许我们通过借用(borrowing)来临时访问数据。

引用与解引用

  • 引用:& 符号用来创建一个引用,允许你使用值但不拥有它。
  • 解引用:* 符号用来访问引用指向的实际数据。
1
2
3
4
5
6
7
8
9
10
fn main() {
let s1 = String::from("hello");
let len = calculate_length(&s1); // 将 s1 的引用传递给函数

println!("The length of '{}' is {}.", s1, len);
}

fn calculate_length(s: &String) -> usize { // s 是一个引用
s.len()
} // 这里 s 离开作用域,但因为它只是引用,所以不会影响原始数据

切片类型

切片让你可以引用集合的一部分,而不需要整个集合的所有权。

1
2
3
4
5
6
7
8
9
fn main() {
let s = String::from("hello world");

let hello = &s[0..5]; // 字符串切片
let world = &s[6..11];

println!("First word: {}", hello);
println!("Second word: {}", world);
}

结论

Rust 的所有权系统是一种独特的方法,用于解决内存管理和并发问题,同时保持性能。虽然一开始可能看起来有些复杂,但它极大地减少了运行时错误,提高了代码的安全性和可靠性。

使用Rust编写的程序由于是编译型的语言,其速度远超Python等脚本程序,所以本文尝试使用Rust编写一段\(f(x)=x^2\)在区间\((0,10)\)上的值,并将其保存在function_data.csv文件中,以供veusz调用,绘制函数图像。

编写Rust程序

  1. 建立项目: drawfunction

    1
    2
    cargo new drawfunction
    cd drawfunction

  2. 确保Cargo.toml文件中添加了必要的依赖项:

    1
    2
    [dependencies]
    csv = "1.1"

  3. 编辑源码:

    drawfunction/src/main.rs
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    use csv::Writer;
    use std::error::Error;

    fn main() -> Result<(), Box<dyn Error>> {
    let mut wtr = Writer::from_path("function_data.csv")?;

    // 写入CSV头部
    wtr.write_record(&["x", "x_squared"])?;

    // 生成数据并写入CSV
    for x in 0..=100 { // 使用0到10之间的101个点
    let x_val = (x as f64) / 10.0;
    let y_val = x_val.powi(2);
    wtr.write_record(&[x_val.to_string(), y_val.to_string()])?;
    }

    wtr.flush()?;
    println!("CSV data has been written.");
    Ok(())
    }

    使用cargo run运行这个程序将会在当前目录下生成一个名为function_data.csv的文件,其中包含了我们计算的数据。

  4. 编译程序:

    1
    2
    cargo build --release
    ./drawfunction/target/release/montecarlo

    第1行为编译程序,第2行为运行这个程序. 可以明显的感觉到程序运行时,瞬间在当前目录下生成一个名为function_data.csv的文件,其中包含了我们计算的数据。

使用Veusz绘图

  1. 启动veusz: 打开Veusz应用程序。
  2. 导入数据:通过“数据->导入”菜单选项导入你用Rust生成的function_data.csv文件。
  3. 创建图表:
  • Veusz中,选择Insert-->add Graph以创建一个新的图形窗口。
  • 在工具栏中选择Plot points with lines and eorrbars.
  • 选中上述组件,在X data 中选择 x , Y data 中选择x_squared.

简介

Rust是一种由Mozilla开发的系统级编程语言,专注于速度、内存安全和并行处理。它通过所有权模型、借用检查器及生命周期概念,在编译时确保内存安全,有效避免了空指针异常、悬挂指针和数据竞争等问题。Rust支持高效的并发编程,并提供了高层次抽象而无运行时开销,保证程序性能。其自带的包管理和构建工具Cargo简化了依赖管理和构建流程。Rust适用于多种场景,包括操作系统、嵌入式系统、WebAssembly、游戏开发及高性能网络服务。凭借现代化的语法、强大的工具链支持(如rustfmt和clippy)以及跨平台兼容性,Rust在追求高效和安全的项目中展现出色。活跃的开源社区和丰富的学习资源进一步增强了它的吸引力,使其成为开发下一代可靠软件的理想选择。

创建Rust项目

  1. 安装Rust: 各种系统中Rust的安装,请参考Rust官网说明. 在ArchLinux环境下,安装命令为

    1
    sudo pacman -S rust

  2. 创建新项目

    1
    2
    cargo new hello_world
    cd hello_world

    运行cargo new hello_world后,会建立一个文件夹hello_world, 其内部包含:Cargo.lock,Cargo.toml,src,target.

  3. 编辑src/main.rs: 打开src/main.rs文件,并添加一些简单的Rust代码,比如打印斐波那契数列的前10项:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    fn fibonacci(n: u32) -> u32 {
    match n {
    0 => 0,
    1 => 1,
    _ => fibonacci(n - 1) + fibonacci(n - 2),
    }
    }

    fn main() {
    for i in 0..10 {
    println!("fibonacci({}) = {}", i, fibonacci(i));
    }
    }

  4. 编译并运行项目:

    1
    cargo run

    运行上述命令后,将输出斐波那契数列的前10项。

使用EVCXR进行交互式探索

为了能够更灵活地试验代码片段而不必每次都通过cargo run执行整个程序,我们可以使用evcxr,这是一个为Rust设计的Jupyter内核,也可以作为独立的REPL使用。

  1. 安装EVCXR:

    1
    cargo install evcxr

  2. 启动EVCXR

    1
    evcxr

  3. EVCXR中试验代码:一旦进入了evcxr环境,您可以直接输入Rust表达式或小段代码进行试验,例如尝试上面定义的斐波那契函数:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    // 定义斐波那契函数
    :dep num_traits

    fn fibonacci(n: u32) -> u32 {
    match n {
    0 => 0,
    1 => 1,
    _ => fibonacci(n - 1) + fibonacci(n - 2),
    }
    }

    // 测试斐波那契函数
    for i in 0..10 {
    println!("fibonacci({}) = {}", i, fibonacci(i));
    }

    请注意,在evcxr中,您可能需要先声明对外部依赖项(如果有的话)使用:dep命令。在这个特定的例子中,我们实际上不需要额外的依赖项,所以可以忽略:dep num_traits这一行。 通过这种方式,您可以在不离开终端的情况下快速测试Rust代码片段,这对于学习语言特性、调试问题或简单地验证想法都非常有用。此外,由于evcxr支持类似Jupyter notebook的交互方式,它也适合用于教学和演示目的。

Rust代码编译

  1. 确保处于正确的目录: 请确保您位于包含Cargo.toml文件的项目根目录下。这个文件定义了项目的配置、依赖和其他相关信息。在上述例子中,就是位于目录hello_world下.

  2. 编译项目:

    1
    cargo build

    这将根据当前环境(调试或发布)编译项目。默认情况下,cargo build会在调试模式下进行编译,生成的二进制文件会位于target/debug/目录下。

  3. 编译优化版本(发布模式):

    编译优化后的版本(适用于生产环境)
    1
    cargo build --release

    这将在target/release/目录下生成经过优化的二进制文件。相比调试模式,发布模式下的编译会花费更多时间,但是生成的二进制文件通常更小且执行速度更快。

  4. 直接运行可执行文件:

  • 对于调试模式编译的二进制文件,可以直接从target/debug/目录下运行:
    1
    ./target/debug/your_project_name
  • 对于发布模式编译的二进制文件,则从target/release/目录下运行:
    1
    ./target/release/your_project_name

    这里的your_project_name是指您Cargo项目的名称,也就是Cargo.toml文件中的[package] name字段所指定的名字。

蒙特卡洛积分是一种通过随机抽样来估计定积分的方法。我们可以用它来近似计算函数 \(y=x^2\)在某个区间上的积分,并基于此结果绘制出其积分函数的图像。 为了验证蒙卡的计算效果,我们在同一张图中给制出蒙卡积分结果和解析解析的结果,对比结果可以证实蒙特卡洛方法的可行性。代码如下:

Comparison between Monte Carlo Integration and Analytical Solutions
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
#! /usr/bin/env python3
# vim:fenc=utf-8
import numpy as np
import matplotlib.pyplot as plt

# 定义蒙特卡洛积分函数
def monte_carlo_integration(func, a, b, n=1000):
x_samples = np.random.uniform(a, b, n)
y_samples = func(x_samples)
# 计算平均高度
average_height = np.mean(y_samples)
# 计算积分值
integral = average_height * (b - a)
return integral

# 定义被积函数
def f(x):
return x**2

# 解析解计算积分函数
def analytic_integral_f(x,a):
return (x**3)/3-(a**3)/3

# 使用蒙特卡洛方法求解积分
a, b = -5, 5 # 积分区间
x_points = np.linspace(a, b, 100) # 用于绘图的x点
monte_carlo_values = []

for i in range(len(x_points)):
# 对每个x点计算从a到该点的积分
integral_value = monte_carlo_integration(f, a, x_points[i], n=20000)
monte_carlo_values.append(integral_value)

# 计算解析解
analytic_values = analytic_integral_f(x_points,a)

# 绘制积分后的函数图像
plt.figure(figsize=(10, 8))
plt.plot(x_points, monte_carlo_values, label='Monte Carlo Integral of $x^2$', color='blue', linestyle='--')
plt.plot(x_points, analytic_values, label='Analytical Integral of $x^2$', color='red')

# 添加标题和标签
plt.title('Comparison between Monte Carlo and Analytical Integration of $x^2$')
plt.xlabel('x')
plt.ylabel('Integral Value')

# 显示图例
plt.legend()

# 显示网格
plt.grid(True)

# 显示图像
plt.show()