我正在try 将matlab中的计算转换为python.这是matlab中的代码:

% Solving for the bandpass correction

T = 8050; % temperature in kelvin
c = 2.99e10; % speed of light in cm/s
h = 6.626e-27; % planck constant in ergs
k = 1.38e-16; % boltzmann constant in erg/K
x1 = 3e-4; % lower wavelength in cm
x2 = 13e-4; % upper wavelength in cm

fun = @(x) ((2 * h * c^2) ./ x.^5) ./ (exp((h * c) ./ (x * k * T)) - 1); % planck function
f_deriv = @(x) ((2 * h^2 * c^3) ./ (x.^6 * k)) .* (exp((h .* c) ./ (x* k * T)) ./ (T * (exp((h * c) ./ (x * k * T)) - 1).^2));
numerator = integral(f_deriv,x1,x2);
denominator = (4 * integral(fun,x1,x2));
beta = numerator / denominator;
fprintf('The numerator is %f \n',numerator)
fprintf('The denominator is %f \n',denominator)
fprintf('The bandpass value is %f \n',beta)

下面是我用python编写代码的try :

# Solving for the bandpass correction:

from scipy.integrate import quad
import numpy as np

T = 8050  # temperature in kelvin
c = 2.99e10  # speed of light in cm/s
h = 6.626e-27  # planck constant in ergs
k = 1.38e-16  # boltzmann constant in erg/K
x1 = 3e-4  # lower wavelength in cm
x2 = 13e-4  # upper wavelength in cm


def p_function(x):
    return ((2 * h * c ** 2) / x ** 5) / (np.exp((h * c)/(x * k * T)) - 1)  # The Planck function


def p_deriv(x):
    return ((2 * h ** 2 * c ** 3) / (x ** 6 * k)) * (np.exp((h * c)/(x * k * T))/(T * (np.exp((h * c) / (x * k * T)) - 1) ** 2))


numerator = quad(p_deriv, x1, x2)
denominator = 4 * quad(p_function, x1, x2)
beta = numerator[0] / denominator[0]
print("The numerator is", numerator[0])
print("The denominator is", denominator[0])
print("The bandpass value is", beta)

两者之间的β值不一致,我看不出是什么导致了差异.我将非常感谢任何帮助来协调这一点.谢谢

推荐答案

通常在翻译MATLAB时,正确的形状/尺寸很重要.但当我用八度音阶运行代码时,我看到所有变量都是(1,1),"标量".所以尺寸不应该是个问题.

让我们判断函数值:

>> fun(1)
ans =  0.066426
>> f_deriv(1)
ans =  0.066432

您的numpy个值看起来相同:

In [3]: p_function(1)
Out[3]: 0.06642589646577152
In [4]: p_deriv(1)
Out[4]: 0.06643181982384715

整合结果:

>> numerator
numerator =  795765635.47589

In [7]: numerator = quad(p_deriv, x1, x2)
In [8]: numerator
Out[8]: (795765635.4758892, 0.030880182071769013)

阅读文档,我们看到quad返回(默认)两个值,一个是整数,一个是误差估计.元组的第一个匹配MATLAB.你用numerator[0].

denominator相同-在乘法之前,注意从元组中提取值:

In [10]: denominator = 4 * quad(p_function, x1, x2)[0]
In [11]: denominator
Out[11]: 2568704689.659281
In [12]: numerator[0] / denominator
Out[12]: 0.309792573151506

>> beta
beta =  0.30979

如果不进行校正,结果将是4倍大

In [13]: denominator = 4 * quad(p_function, x1, x2)
In [14]: denominator
Out[14]: 
(642176172.4148202,
 0.006319781838840299,
 642176172.4148202,
 0.006319781838840299,
 642176172.4148202,
 0.006319781838840299,
 642176172.4148202,
 0.006319781838840299)
In [15]: numerator[0] / denominator[0]
Out[15]: 1.239170292606024

有时(甚至经常)我们需要一步一步地测试值.即使从头开始开发numpy个代码,我也会测试所有步骤.

Python相关问答推荐

为什么我的(工作)代码(生成交互式情节)在将其放入函数中时不再工作?

使用imap-tools时错误,其邮箱地址包含域名中的非默认字符

删除pandas rame时间序列列中未更改的值

使用Beautiful Soup获取第二个srcset属性

如何修复使用turtle和tkinter制作的绘画应用程序的撤销功能

将轨迹优化问题描述为NLP.如何用Gekko解决这个问题?当前面临异常:@错误:最大方程长度错误

Python daskValue错误:无法识别的区块管理器dask -必须是以下之一:[]

滚动和,句号来自Pandas列

ModuleNotFound错误:没有名为Crypto Windows 11、Python 3.11.6的模块

在Google Colab中设置Llama-2出现问题-加载判断点碎片时Cell-run失败

如何在Django基于类的视图中有效地使用UTE和RST HTIP方法?

如何在给定的条件下使numpy数组的计算速度最快?

使用Python和文件进行模糊输出

在matplotlib中删除子图之间的间隙_mosaic

如何使用Numpy. stracards重新编写滚动和?

寻找Regex模式返回与我当前函数类似的结果

循环浏览每个客户记录,以获取他们来自的第一个/最后一个渠道

在Python中从嵌套的for循环中获取插值

用两个字符串构建回文

为什么后跟inplace方法的`.rename(Columns={';b';:';b';},Copy=False)`没有更新原始数据帧?