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}")

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}")

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() 

# 导入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}") 

#! /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}")

#! /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")

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

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

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

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

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}")