MENU

计算级数

• March 21, 2024 • Read: 1394 • Partial Differential Equation,Python,Numerical Analysis

PDE采用解析方法求解后,通常解的形式都是傅里叶级数形式表示的,计算各个采样点的值需要完成级数的求和计算。

$$ u(z,t) = \frac{4 u_0}{\pi} \sum_{n=1}^{\infty} \frac{1}{2n-1} \sin(\frac{2n-1}{2H} \pi z) \exp(- \frac{((2n-1)\pi}{2H})^2 C_v t) $$

一维空间中的非稳态问题,场函数必定是一个二元函数,自变量有$x$和$t$,级数的计算还需要考虑$n$的取值,需要计算每一个$n$对应下的结果,因此,需要不断地更新$n$生成新的函数,然后求和就可以得到$z,t$对你应下的结果。

利用sympy实现计算,其中summation函数能够实现级数的计算,因此,只需要每一次更新$z,t$的值即可。

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt


n = sp.symbols('n')
H = 1.0
zv = np.linspace(0, H, 11)
t = 1e5
tv = np.linspace(0, t, 50)
cv = 1e-5
u0 = 1000


# 初始化结果矩阵,时间为竖向,空间值按照列进行放置
sol = np.zeros([tv.shape[0], tv.shape[0]])

for row in range(tv.shape[0]):
    for col in range(zv.shape[0]):
        f = 1 / (2 * n - 1) * sp.sin((2 * n - 1)/ (2 * H) * sp.pi * zv[col]) * sp.exp(-1 * ((2*n-1)*sp.pi/ (2*H))**2 * 1e-5 * tv[row])
        sol[row, col] = 4 * u0 / np.pi * sp.sumation(f, (n, 1, 50)).evalf()

evalf()方法,能够将计算结果转换为浮点数进行表示,sympy中类似的还有sympify和lambdify,都能实现将结果以数值的形式表示出来。

Archives Tip
QR Code for this page
Tipping QR Code