我试图计算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
)要么
你的代码差不多就是这样
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的intégral函数可以完成它的工作。
一般来说,接受的答案是我所说的最好的方法,但如果你的功能有某些限制是允许的,那么还有第二种方法。
有两个函数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));