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

Evaluate indefinite integral numerically in matlab/mathematica that it cannot do symbolically

我正在尝试在 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 循环中的两个计算都给出了 fa 的粗略近似值,而不是它们乘以 dt 后的积分。

在数学堆栈交换中,我发现了一个 question 可以推导出一个点的积分的有限差分方法。然而,当我在 Matlab 中进行计算时,它输出了 f 的缩小版本,这在绘图后很明显(请参阅上文了解我所说的缩小的意思)。我认为这是因为对于较小的间隔,积分基本上以不同的精度逼近函数(再次参见上文)。

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

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

求积分
  1. 除了默认的积分计算器,我可以直接用积分计算器求积分吗? (int,integral,trapz)

  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 积分函数可以正常工作。

例如让我们拿 = the function above will compute if you wish to compute 把上面的 0 替换成你喜欢的常量 a

现在 and so you should get F containing a discretization of

到目前为止,公认的答案是我认为最好的方法,但如果允许对您的功能进行某些限制,那么还有第二种方法。

两个函数 fg 见下文

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));