how does double integration work
显示 更早的评论
Hello eveyone, please I am trying to find the first and second integral of a given fuction (say x). I wrote a code below. But it seems the second integration plot does not work appropriately. pls can someone help me. And also make other neccessary corrections. Thank you for your help
x = -100: 1/10:100;
f = 6*x;
first_integral_plot = cumtrapz(x, f); % Integration
second_integral_plot = cumtrapz(x, first_integral_plot ); % Integration
figure;subplot; subplot(311);plot(real(f));title('original function')
subplot(312);plot(first_integral_plot);title('integration2')
subplot(313); plot(second_integral_plot);title('integration2')
采纳的回答
x = -100: 1/10:100;
f = 6*x;
first_integral_plot = cumtrapz(x, f); % Integration
second_integral_plot = cumtrapz(x, first_integral_plot ); % Integration
Why do you think this has not worked? The first integral is of a linear kernel. So we should see quadratic behavior in the output.
plot(x,first_integral_plot)
That looks nicely quadratic to me. Now, the integral of a quadratic would seem to be cubic in nature.
plot(x,second_integral_plot)
I'l bite. What seems wrong to you?
Let me guess. You forget that an integral actually has an undetermined constant in there. So the indefinite integral of 6*x should be 3*x^2 + C. Of course, cumtrapz, being a purely numerical code, does not understand what to do with that undetermined constant. So it sort of ignores it. This is NOT the same thing as to asume the constant was zero. In fact, the implicit presumption of cumtrapz is that at the first point, the integral was ZERO. And that is correct, because cumtrapz is a DEFINITE integration tool. That is, what is the integral of F from -100 to some point X? When X is -100, that value is ZERO.
So what is cumtrapz producing? I'll do the same computation in a symbolic form, and then plot what we get.
syms x
F = 6*x;
Fint = int(F,x,-100,x)
Again, that is a DEFINITE INTEGRAL. Does it look like what cumtrapz produced?
fplot(Fint,[-100,100])
Now, when you integrate that a second time, we see:
Fint2 = expand(int(Fint,x,-100,x))
fplot(Fint2,[-100,100])
It should not be surpriise that the two agree remarkably well. And of course, we can check the result:
diff(Fint2,x,2)
All seems to be happy.
Never forget those constants of integration. Sometimes they come back and bite you.
2 个评论
Eshovo Moses
2021-3-5
I'm confused. Why would I want to write code to do what I already know works? NEVER WRITE CODE TO DO WHAT ALREADY EXISTS IN THE FORM OF HIGH QUALITY CODE WRITTEN BY PROFESSIONAL PROGRAMMERS.
In this case, I've already shown two fundamental ways to accomplish what you asked about, using existing tools. They are perfectly consistent in their results. What matters here is to understand the mathematics, and I explained that. In fact, there are many ways to solve the problem you are asking to solve, all already available.
So if you want a third solution, consider what it means to solve the ordinary differential equation
y''(x) = F(x) = 6*x
over the interval [-100,100]. Supply a pair of initial conditions for y(-100), and y'(-100). (Hint: you might decide to use all sorts of methods to solve that ODE, from Euler's method to a Runge-Kutta code. But far better is to follow the advice I gave above. Use the existing solution codes in MATLAB. The tools are already there. You might start with ODE45. dsolve would also apply.)
Lets see, does this produce the same solution as I found above?
syms x y(x)
F = 6*x;
dy(x) = diff(y,x);
ysol = dsolve(diff(y,x,2) == F,y(-100) == 0,dy(-100) == 0)
Of course it does. But my point is, no, I won't write an ode solving code for you. You already have many of them, certainly written better than I would in a limited amount of time.
If you want some code, then feel free to write it. If you don't understand the coding involved, then attempting to write that code is a good way to learn. Go for it. Your choice.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Numeric Solvers 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
