关于#python#的问题：求一个二维和三维的对流扩散方程Python的程序 要求有误差分析和解析解和数值解的比较 事成六张

求一个二维和三维的对流扩散方程Python的程序 要求有误差分析和解析解和数值解的比较 事成六张

对流扩散方程是一类常见的偏微分方程，其数值求解方法有很多种，包括有限差分方法、有限元方法、谱方法等。下面是一个使用有限差分方法求解二维和三维对流扩散方程的 Python 程序，并包含误差分析和解析解和数值解的比较。

二维对流扩散方程的数值求解：

import numpy as np
import matplotlib.pyplot as plt

# 空间区间和步长设置
Lx, Ly = 1.0, 1.0
dx, dy = 0.01, 0.01
Nx, Ny = int(Lx/dx), int(Ly/dy)

# 时间区间和步长设置
T, dt = 1.0, 0.001
Nt = int(T/dt)

# 物理参数设置
rho, mu, kappa = 1.0, 0.001, 0.001

# 初始条件设置
u0 = np.zeros((Nx, Ny))
u0[int(Nx/4):int(3*Nx/4), int(Ny/4):int(3*Ny/4)] = 1.0

# 数值解计算
u = u0.copy()
for n in range(Nt):
    u[1:Nx-1,1:Ny-1] += dt*(-rho*u[1:Nx-1,1:Ny-1]*(u[2:Nx,1:Ny-1]-u[0:Nx-2,1:Ny-1])/dx
                            -rho*u[1:Nx-1,1:Ny-1]*(u[1:Nx-1,2:Ny]-u[1:Nx-1,0:Ny-2])/dy
                            +mu*(u[2:Nx,1:Ny-1]-2*u[1:Nx-1,1:Ny-1]+u[0:Nx-2,1:Ny-1])/dx**2
                            +mu*(u[1:Nx-1,2:Ny]-2*u[1:Nx-1,1:Ny-1]+u[1:Nx-1,0:Ny-2])/dy**2)

# 解析解计算
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
X, Y = np.meshgrid(x, y)
u_exact = np.exp(-2*rho/kappa*T)*np.sin(np.pi*X)*np.sin(np.pi*Y)

# 误差分析和比较
error = np.linalg.norm(u_exact - u, ord=2)/(Nx*Ny)
print('L2误差：', error)
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.imshow(u_exact, cmap='jet', extent=[0, Lx, 0, Ly])
plt.colorbar()
plt.title('解析解')
plt.subplot(1, 2, 2)
plt.imshow(u, cmap='jet', extent=[0, Lx, 0, Ly])
plt.colorbar()
plt.title('数值解')
plt.show()

三维对流扩散方程的数值求解：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 空间区间和步长设置
Lx, Ly, Lz = 1.0, 1.0, 1.0
dx, dy, dz = 0.01, 0.01, 0.01
Nx, Ny, Nz = int(Lx/dx), int(Ly/dy), int(Lz/dz)

# 时间区间和步长设置
T, dt = 1.0, 0.001
Nt = int(T/dt)

# 物理参数设置
rho, mu, kappa = 1.0, 0.001, 0.001

# 初始条件设置
u0 = np.zeros((Nx, Ny, Nz))
u0[int(Nx/4):int(3*Nx/4), int(Ny/4):int(3*Ny/4), int(Nz/4):int(3*Nz/4)] = 1.0

# 数值解计算
u = u0.copy()
for n in range(Nt):
    u[1:Nx-1,1:Ny-1,1:Nz-1] += dt*(-rho*u[1:Nx-1,1:Ny-1,1:Nz-1]*(u[2:Nx,1:Ny-1,1:Nz-1]-u[0:Nx-2,1:Ny-1,1:Nz-1])/dx
                            -rho*u[1:Nx-1,1:Ny-1,1:Nz-1]*(u[1:Nx-1,2:Ny,1:Nz-1]-u[1:Nx-1,0:Ny-2,1:Nz-1])/dy
                            -rho*u[1:Nx-1,1:Ny-1,1:Nz-1]*(u[1:Nx-1,1:Ny-1,2:Nz]-u[1:Nx-1,1:Ny-1,0:Nz-2])/dz
                            +mu*(u[2:Nx,1:Ny-1,1:Nz-1]-2*u[1:Nx-1,1:Ny-1,1:Nz-1]+u[0:Nx-2,1:Ny-1,1:Nz-1])/dx**2
                            +mu*(u[1:Nx-1,2:Ny,1:Nz-1]-2*u[1:Nx-1,1:Ny-1,1:Nz-1]+u[1:Nx-1,0:Ny-2,1:Nz-1])/dy**2
                            +mu*(u[1:Nx-1,1:Ny-1,2:Nz]-2*u[1:Nx-1,1:Ny-1,1:Nz-1]+u[1:Nx-1,1:Ny-1,0:Nz-2])/dz**2)

# 解析解计算
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
z = np.linspace(0, Lz, Nz)
X, Y, Z = np.meshgrid(x, y, z)
u_exact = np.exp(-3*rho/kappa*T)*np.sin(np.pi*X)*np.sin(np.pi*Y)*np.sin(np.pi*Z)

# 误差分析和比较
error = np.linalg.norm(u_exact - u, ord=2)/(Nx*Ny*Nz)
print('L2误差：', error)
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.plot_surface(X, Y, u_exact[:,:,int(Nz/2)], cmap='jet')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u')
ax.set_title('解析解')
ax = fig.add_subplot(1, 2, 2, projection='3d')
ax.plot_surface(X, Y, u[:,:,int(Nz/2)], cmap='jet')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u')
ax.set_title('数值解')
plt.show()

这里使用了 NumPy 和 Matplotlib 库来进行数值计算和可视化，其中 np.linalg.norm 函数用于计算 L2 范数，plt.imshow 函数用于显示二维图像，Axes3D.plot_surface 函数用于显示三维图像。通过比较解析解和数值解的可视化效果和 L2 误差，可以评估数值求解的精度和准确性。

  • 这有个类似的问题, 你可以参考下: https://ask.csdn.net/questions/7660763
  • 你也可以参考下这篇文章：如何用python实现高斯分布
  • 除此之外, 这篇博客: 机器学习之一元线性回归模型梯度下降算法的Python实现中的 机器学习之一元线性回归模型梯度下降算法的Python实现 部分也许能够解决你的问题, 你可以仔细阅读以下内容或跳转源博客中阅读:

  • 前言： 关于一元线性回归梯度下降法的详细理论知识可至我机器学习类别的博文中查看，本文基于Python实现梯度下降算法。

    一、算法实现代码
    本数据的学习速率选择0.0001，迭代次数为50次，初始参数值为0

    import numpy as np
import matplotlib.pyplot as plt

# 导入数据
data = np.genfromtxt('data.csv', delimiter=',')
x_data = data[:, 0]
y_data = data[:, 1]

# 学习率
lr = 0.0001
# 一元线性模型: y_data = theta0 + theta1 * x_data
theta0_init = 0
theta1_init = 0
# 梯度下降法最大迭代次数
ite = 50

# 代价函数
def compute_error(theta0, theta1, x_data, y_data):
    totleError = 0
    for i in range(0, len(x_data)):
        totleError += ((theta0 + theta1 * x_data[i]) - y_data[i]) **2
    return totleError / float(len(x_data)) / 2.0


# 梯度下降
def gradient_descent_runner(theta0, theta1, lr, ite, x_data, y_data):
    m = float(len(x_data))  # 计算样本数量

    functionFigure_num = 1  # 图形个数
    # 循环迭代进行梯度下降
    for i in range(ite):
        theta0_grad = 0
        theta1_grad = 0
        for j in range(0, len(x_data)):
            theta0_grad += 1 / m * ((theta0 + (theta1 * x_data[j])) - y_data[j])
            theta1_grad += 1 / m * ((theta0 + (theta1 * x_data[j])) - y_data[j]) * x_data[j]

        # 同步更新theta0 和theta1
        theta0 = theta0 - (lr * theta0_grad)
        theta1 = theta1 - (lr * theta1_grad)

        # 每迭代5次显示一次当前梯度下降的一元线性模型图像
        if i % 5 == 0:
            # 在一个窗口中每行5个显示图像
            print('iteration:', i)
            plt.subplot(2, 5, functionFigure_num)
            plt.plot(x_data, y_data, 'b.')
            line, = plt.plot(x_data, theta0 + (theta1 * x_data), 'r')
            plt.legend(handles=[line], labels=['iterations' + str(functionFigure_num * 5)], loc='best')
            functionFigure_num = functionFigure_num + 1
    plt.show()
    return theta0, theta1

if __name__ == '__main__':
    print('Starting theta0 = {0}, theta1 = {1}, error = {2}'.format(theta0_init, theta1_init, compute_error(theta0_init, theta1_init, x_data, y_data)))
    print('Gradient Descent Running...')
    theta0_end, theta1_end = gradient_descent_runner(theta0_init, theta1_init, lr, ite, x_data, y_data)
    print('After {0} iterations theta0 = {1}, theta1 = {2}, error = {3}'.format(ite, theta0_end, theta1_end, compute_error(theta0_end, theta1_end, x_data, y_data)))

    # 绘制最终梯度下降后的直线
    # plt.plot(x_data, y_data, 'b.')
    # plt.plot(x_data, theta0_end + theta1_end * x_data, 'r')
    # plt.show()

    二、实现结果

    在这里插入图片描述
    由上图可知迭代到最后的线性模型较好的拟合了样本数据

    三、数据下载地址
    链接：https://pan.baidu.com/s/1KhPIzejxFZfbkIGZo8J8lg
    提取码：slg1

  • 您还可以看一下 李飞老师的Python数据分析与爬虫课程中的 Python高级特性：迭代器与生成器小节, 巩固相关知识点