Find integral upper limit knowing the result and having an array
3 次查看(过去 30 天)
显示 更早的评论
Hello,
I am trying to find the upper limit of an integral, knowing the results, but I have an array instead of a function.
If I had a function, for example -2t, and an integral result of lets say -4, i could find the upper limit as shown below, using fzero.
fun = @(t)-2*t;
intfun = @(x)integral(fun,0,x) + 4;
xval = fzero(intfun,[0,10]);
and i find an xval=2 for this case.
Howover now I do not have a function, but an array, for exaple, y2 as shown below. I have different values for y2 as I find y2 solving a single degree of freedom system and it represents the velocity values for each time instance. I need to find the time tout (upper limit of integral) and lets say the results is 10. I tried to solve it as in the example above, however i get an error.
y2=0:0.1:15;
fun_1= y2;
vale1= @(tout)integral(fun_1,0,tout,'ArrayValued',true)+10;
ans2=fzero(vale1,[0,10]);
here is the error:
Error using fzero>localFirstFcnEval (line 729)
FZERO cannot continue because user-supplied function_handle ==> @(tout)integral(fun_1,0,tout,'ArrayValued',true)+10
failed with the error below.
First input argument must be a function handle.
Could you help me with it please?
Thank you!
Angie
2 个评论
回答(2 个)
John D'Errico
2023-1-24
A basic rule is, you CANNOT use integral OR fzero with data, thus a vector of numbers.
You CAN approximate your data using a function, perhaps a spline interpolant, and then integrate that, and use fzero on the result. There is a fundamental difference here that you need to understand.
A second approach would be to form a cumulative integral, perhaps using cumtrapz. And then interpolate that.
For example, can you find the solution of int(sin(x)), on limits of [0,u], such that the integral is 1.75? So you would want to solve for u.
The problem is, I'll give you only data.
x = linspace(0,pi);
y = sin(x);
How would we solve it? A simple solution is as I said, to use cumtrapz.
yint = cumtrapz(x,y);
inttarget = 1.75;
plot(x,yint,'b-')
yline(inttarget,'r')
I''ve plotted the cumulative integral there, on the interval [0,pi].
The left endpoint of the integral is fixed at 0. The right endpoint is unknown. So all we need to do now is find the value of u, such that yint is exactly 1.5. Of course we will never be lucky enough that one of those data points will be exactly 1.75.
Here, since the function is monotone increasing on the interval, we can just swap x and y, and then use interp1.
format long g
u = interp1(yint,x,inttarget,'spline')
How well did we do? We can check against the analytical result.
syms X U
Usol = solve(int(sin(X),X,[0,U]) == inttarget)
vpa(Usol)
The second solution is the one we want. As you can see, the cumtrapz solution was not terrible. We got almost 4 correct digits in the solution.
The other approach may be a little more accurate, since it does not use a trapezoidal integration.
spl = spline(x,y); % Interpolate the curve first with a spline
splint = fnint(spl); % integrate
% Now solve for the upper limit of integration
[u,fval] = fzero(@(u) fnval(splint,u)-fnval(splint,0)-inttarget,[0,pi])
This approach yields a few more significant digits. I would expect that to be the case as long as the data is not noisy. With noisy data, I would use other approaches that would probably involve smoothing the data first.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Discrete Math 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!