在 matlab/mathematica 中以数值方式评估不定积分，它不能以符号方式进行答案

【问题标题】：Evaluate indefinite integral numerically in matlab/mathematica that it cannot do symbolically在 matlab/mathematica 中以数值方式评估不定积分，它不能以符号方式进行
【发布时间】：2018-09-19 03:28:43
【问题描述】：

我正在尝试在 Matlab 和 Mathematica 中计算软件无法以符号方式执行的函数的积分。

到目前为止，这是我的 MatLab 代码，但我知道它可能不是很有帮助。

f = @(t) asin(0.5*sin(t));
a = @(t) sin(t);
F = int(f,t)   % Matlab can't do this
F = 
int(asin(sin(t)/2), t)
A = int(a,t)   % This works
A =
-cos(t)

dt = 1/(N-1); % some small number
for i=1:N
    F(i) = integral(f,(i-1)*dt,i*dt);
    A(i) = integral(a,(i-1)*dt,i*dt);
end

for 循环中的两个计算都给出了f 或a 的粗略近似值，而不是它们乘以dt 后的积分。

在数学堆栈交换中，我发现了一个question，它为一点积分推导了一个有限差分之类的方法。但是，当我在 Matlab 中进行计算时，它会输出一个缩小版本的f，这在绘制后很明显（见上文，我所说的缩小是什么意思）。我认为这是因为对于较小的间隔，积分基本上以不同程度的准确度近似函数（再次参见上文）。

我正在尝试获得积分的符号方程，或者每个位置的函数积分的近似值。

所以我的问题是如果我有一个函数 f ，MatLab 和 Mathematica 不能轻易地取积分

除了默认的积分计算器，我可以直接用积分计算器近似积分吗？ (int,integral,trapz)

或

我可以先用有限差分逼近函数，然后符号地求积分吗？

【问题讨论】：

你是什么意思“除了默认的” - 你将什么归类为默认求解器，为什么你认为它们不起作用？最简单的数值积分器是 Matlab 中的trapz，为什么这不起作用？如果您的限制是不确定的（在这种情况下您不能用数字评估任何东西），您的预期输出是什么？
我解释了为什么它们不起作用，积分间隔太小，除了您输入计算器的函数之外，无法输出任何东西。
请附上minimal reproducible example。我们无法运行您当前的代码来产生您得到的任何非结果，并且您没有建议预期的输出可能是什么。你能举一个例子，其中被积函数可以代数积分，所以我们可以验证任何结果？
我不知道为什么你不能运行提供的代码，这并不难。但我做了一些修改，应该会让你更容易。

标签： matlab wolfram-mathematica

【解决方案1】：

您的代码几乎没问题，只是

for i=1:N
    F(i) = integral(f,0,i*dt);
end

你也可以

F(1)=integral(f,0,dt)
for i=2:N
    F(i) = F(i-1)+integral(f,(i-1)*dt,i*dt);
end

第二种选择肯定更有效

因为原语实际上是 F(x)=int(f(x), 0, x) （0 定义了某个常数）并且对于足够小的 dx，您已经证明 f(x)=int(f(x ), x,x+dx)/dx i。您已经证明 MATLAB 积分函数可以发挥作用。

例如，假设= ，上面的函数将计算，如果你想计算，只需将上面的0替换为你喜欢的常量a。

现在是，所以你应该得到包含离散化的F

【讨论】：

我知道int“它的工作”。也许我的方法有问题，但如果我整合 f(x) 积分应该像转换一样。该函数不是指数函数，因此 f(x) 不等于 int(f,x)。
您是否尝试过修复我拥有您（将 (i-1)dt 替换为 0）？我只是说当 dx 变为零时 f(x)=int(f(x), x,x+dx)/dx 这是你在没有修复的情况下观察到的。
不，我没有，但原因是因为我想移动积分的间隔。从0:dt 开始，直到T-dt:T。有意义吗？
我稍后再试试，谢谢。我暂时不能。
感谢您最终接受这个答案，我明天将删除我的 cmets。

【解决方案2】：

总的来说，公认的答案是我想说的最好的方法，但如果允许对你的功能进行某些限制，那么还有第二种方法。

f 和 g 两个函数见下文

T = 1;  % Period
NT = 1;  % Number of periods
dt = 0.01; % time interval
time = 0:dt:NT*T;  % time

syms t
x = K*sin(2*pi*t+B);   % edit as appropriate

% f = A/tanh(K)*tanh(K*sin(2*pi*t+p))
% g = A/asin(K)*asin(K*sin(2*pi*t+p))

找到公式here

f = A1/tanh(K1)*(2^(2*1)-1)*2^(2*1)*bernoulli(2*1)/factorial(2*1)*x^(2*1-1);
% |K1|<pi/2
g = A2/asin(K2)*factorial(2*0)/(2^(2*0)*factorial(0)^2*(2*0+1))*x^(2*0+1);
% |K2|<1

接受的答案没有这样的限制

N = 60;
for k=2:N
    a1 = (2^(2*k)-1)*2^(2*k)*bernoulli(2*k)/factorial(2*k);
    f = f + A1/tanh(K1)*a1*x^(2*k-1);

    a2 = factorial(2*k)/(2^(2*k)*factorial(k)^2*(2*k+1));
    g = g + A2/asin(K2)*a*x^(2*k+1);
end

MATLAB 可以计算 sin^n(t)，因为 n 是一个整数。

F = int(f,t);
phi = double(subs(F,t,time));

G = int(g,t);
psi = double(subs(G,t,time));

【讨论】：