numpy 如何绘制取决于两个变量的SymPy函数？

qcbq4gxm 于 2023-03-18 发布在其他

关注(0)|答案(1)|浏览(125)

因此，我尝试实现下面的函数Function，它应该会产生下面的图：w，l，δ通常是常量，但是函数C_EF以w/δ作为参数，所以我使用SymPy实现了w/δ。
为了科普从0到Inf的求和，我使用了while循环来检查步骤之间是否有足够小的增量以停止求和。然而，在我的整体实现中，我无法得到与所包含的图形相同的图形。我不确定我的代码哪里出错了，但我正在考虑，依赖于δ^2和w的函数也可能导致w/δ旁边的问题。编辑：下面的屏幕截图显示了我的代码生成的图：My Plot

### IMPORTS AND ABBREVATIONS ###
import math
import numpy as np
import scipy
from sympy import *  # plot, re, im, symbols, I
import sympy as sym

mu_r = 1.05  # Relative permeability
mu_0 = 4 * pi * 10 ** -7  # magnetic permeability of free space
mu = mu_0 * mu_r  # H/m; Absolute permeability
sigma = 7.1 * 10 ** 5  # S/m; Electrical conductivity für NdFeB Magnet
f = 2500  # Hz; carrier Frequency
B_h_s = 1.5  # T; peak Flux density

dAG = 2.5 * 10 ** -3  
w_magnet = 60 * 10 ** -3  
h = 5 * 10 ** -3  
l = 3*w_magnet  
pi = scipy.pi

wd = Symbol('wd')

def P_calc(w):
    dlt = (sym.sqrt(1 / (pi * mu * sigma * f)) * sym.sqrt((h + dAG) / h)).evalf()

    # (11) Sums
    def nth_term(a, n_sum):
        lmd_n = (2 * n_sum + 1) * pi / w 
        beta_n = sym.sqrt(lmd_n ** 2 + 2 * I / (dlt ** 2)) 
        if a == 1:
            return ((lmd_n ** 2 - 2 * im(beta_n) ** 2) * re(beta_n) * lmd_n ** 3 * sym.sinh(re(beta_n) * l)) / \
                   ((2 * n_sum + 1) ** 5 * abs(beta_n) ** 6 * (sym.cosh(re(beta_n) * l) + sym.cos(im(beta_n) * l)))
        elif a == 2:
            return ((lmd_n ** 2 + 2 * re(beta_n) ** 2) * im(beta_n) * lmd_n ** 3 * sym.sin(im(beta_n) * l)) / \
                   ((2 * n_sum + 1) ** 5 * abs(beta_n) ** 6 * (sym.cosh(re(beta_n) * l) + sym.cos(im(beta_n) * l)))

    tolerance = 1e-6
    max_iterations = 1000
    sum_1_C_EF = 0
    n = 0
    while True:
        term_1 = nth_term(1, n)
        sum_1_C_EF += term_1
        n += 1
        # Check convergence criteria
        if math.fabs(term_1) < tolerance or n >= max_iterations:
            break

    sum_2_C_EF = 0
    n = 0
    while True:
        term_2 = nth_term(2, n)
        sum_2_C_EF += term_2
        n += 1
        # Check convergence criteria
        if math.fabs(term_2) < tolerance or n >= max_iterations:
            break

    sum_C_EF = sum_1_C_EF + sum_2_C_EF
    C_EF_sinh_cosh = ((sym.cosh(wd) + sym.cos(wd)) / (sym.sinh(wd) - sym.sin(wd)))
    C_EF_term = ((32 * w) / ((pi ** 5) * l)) * (wd ** 3)

    # (11) complete
    C_EF = 1 - C_EF_term * C_EF_sinh_cosh * sum_C_EF

    p = plot(C_EF, (wd, 0, 6), ylabel='C_EF', line_color='b', show=False)
    p.show()

    return P_e

P_e = P_calc(w_magnet)

numpy

来源：https://stackoverflow.com/questions/75715245/how-do-i-plot-a-sympy-function-which-is-depending-on-two-variables

1条答案

按热度按时间

2ul0zpep1#

对于这个任务，我不会使用SymPy，因为re和im函数将应用简化，这实际上使表达式非常长，因此难以计算。
下面是我设置代码的方法：

import numpy as np
import matplotlib.pyplot as plt

mu_r = 1.05  # Relative permeability
mu_0 = 4 * np.pi * 10 ** -7  # magnetic permeability of free space
mu = mu_0 * mu_r  # H/m; Absolute permeability
sigma = 7.1 * 10 ** 5  # S/m; Electrical conductivity für NdFeB Magnet
f = 2500  # Hz; carrier Frequency
B_h_s = 1.5  # T; peak Flux density

dAG = 2.5 * 10 ** -3  
h = 5 * 10 ** -3 
w_magnet = 60 * 10 ** -3  
delta = np.sqrt(1 / (np.pi * mu * sigma * f)) * np.sqrt((h + dAG) / h)
w = w_magnet

# break down pieces for easy of coding and reading
lam = lambda n: (2 * n + 1) * np.pi / w
beta = lambda n: np.sqrt(lam(n)**2 + 2j / delta**2)
beta_r = lambda n: np.real(beta(n))
beta_i = lambda n: np.imag(beta(n))
den = lambda l, n: (2 * n + 1)**5 * np.absolute(beta(n))**6 * (np.cosh(np.float128(beta_r(n) * l)) + np.cos(beta_i(n) * l))

def s1(l, N):
    s = lambda n: (lam(n)**2 - 2 * beta_i(n)**2) * beta_r(n) * lam(n)**3 * np.sinh(np.float128(beta_r(n) * l)) / den(l, n)
    return sum(s(n) for n in range(N + 1))

def s2(l, N):
    s = lambda n: (lam(n)**2 + 2 * beta_r(n)**2) * beta_i(n) * lam(n)**3 * np.sin(beta_i(n) * l) / den(l, n)
    return sum(s(n) for n in range(N + 1))

# sum up to N
N = 100
t1 = lambda l: (32 * w / (np.pi**5 * l))
t2 = lambda w_d: (np.cosh(w_d) + np.cos(w_d)) / (np.sinh(w_d) - np.sin(w_d))
t3 = lambda l: s1(l, N) + s2(l, N)
C_EF = lambda w_d, l: 1 - w_d**3 * t1(l) * t2(w_d) * t3(l)

C_EF = np.vectorize(C_EF)
w_d = np.linspace(0, 6, 10)

plt.figure()
for m in [1, 1.5, 3, 8]:
    plt.plot(w_d, C_EF(w_d, m * w), label=r"$l = %s \cdot w$" % m)
plt.legend()

很明显，这不是您要找的图表。我已经用键入的内容仔细检查了公式：在我看来它们看起来是一样的。
注：
1.我没有求出无穷大，而是求出N=100，这由你来决定，我也试过N=5，结果似乎没有太大变化。
1.我将np.sinh和np.cosh的参数转换为np.float128，因为对于相对较大的n，这些项的结果将不适合常规float64。

赞(0）回复(0）举报 2023-03-18

我来回答

numpy 如何绘制取决于两个变量的SymPy函数？

1条答案

相关问题

热门标签

最新问答