通过MATLAB中的梯形和进行积分 [英] Integration via trapezoidal sums in MATLAB

问题描述

我需要使用梯形和求函数积分的帮助.

I need help finding an integral of a function using trapezoidal sums.

程序应使用n = 1, 2, 3, ...连续获取梯形和 直到出现两个相邻的n值相差小于给定公差的时间间隔.我希望在WHILE循环中至少有一个FOR循环，并且我不想使用trapz函数.该程序有四个输入:

The program should take successive trapezoidal sums with n = 1, 2, 3, ... subintervals until there are two neighouring values of n that differ by less than a given tolerance. I want at least one FOR loop within a WHILE loop and I don't want to use the trapz function. The program takes four inputs:

• f:x函数的函数句柄.
• a:一个实数.
• b:大于a的实数.
• tolerance:一个正数，很小的实数

• f: A function handle for a function of x.
• a: A real number.
• b: A real number larger than a.
• tolerance: A real number that is positive and very small

我遇到的问题是尝试实现梯形和的公式 Δx/2[y0 + 2y1 + 2y2 + … + 2yn-1 + yn]

The problem I have is trying to implement the formula for trapezoidal sums which is Δx/2[y0 + 2y1 + 2y2 + … + 2yn-1 + yn]

这是我的代码，卡住的区域是FOR循环内的求和"部分.我正在尝试总结2y2 + 2y3....2yn-1，因为我已经解释了2y1.我得到了一个答案，但是并没有达到应有的准确度.例如，我得到的是6.071717974723753而不是6.101605982576467.

Here is my code, and the area I'm stuck in is the "sum" part within the FOR loop. I'm trying to sum up 2y2 + 2y3....2yn-1 since I already accounted for 2y1. I get an answer, but it isn't as accurate as it should be. For example, I get 6.071717974723753 instead of 6.101605982576467.

感谢您的帮助！

function t=trapintegral(f,a,b,tol)
format compact; format long;
syms x;
oldtrap = ((b-a)/2)*(f(a)+f(b));
n = 2;
h = (b-a)/n;
newtrap = (h/2)*(f(a)+(2*f(a+h))+f(b));
while (abs(newtrap-oldtrap)>=tol)
    oldtrap = newtrap;
    for i=[3:n]
        dx = (b-a)/n;
        trapezoidsum = (dx/2)*(f(x) + (2*sum(f(a+(3:n-1))))+f(b));
        newtrap = trapezoidsum;
    end
end
t = newtrap;
end

推荐答案

此代码不起作用的原因是，梯形规则的求和中存在两个小错误.我要确切地指的是这样的声明:

The reason why this code isn't working is because there are two slight errors in your summation for the trapezoidal rule. What I am precisely referring to is this statement:

trapezoidsum = (dx/2)*(f(x) + (2*sum(f(a+(3:n-1))))+f(b));

调用梯形积分法则的方程:

Recall the equation for the trapezoidal integration rule:

^{来源:维基百科}

对于第一个错误，在包含起点时，f(x)应该为f(a)，并且不应保留为符号.实际上，您应该简单地删除syms x语句，因为它在脚本中没有用.通过参考上述公式，a对应于x1.

For the first error, f(x) should be f(a) as you are including the starting point, and shouldn't be left as symbolic. In fact, you should simply get rid of the syms x statement as it is not useful in your script. a corresponds to x1 by consulting the above equation.

下一个错误是第二项.实际上，您实际上需要将索引值(3:n-1)乘以dx.另外，这实际上应该来自(1:n-1)，我将在后面解释.上面的公式从2变为N，但是出于我们的目的，当您将代码设置为这样时，我们将从1变为N-1.

The next error is the second term. You actually need to multiply your index values (3:n-1) by dx. Also, this should actually go from (1:n-1) and I'll explain later. The equation above goes from 2 to N, but for our purposes, we are going to go from 1 to N-1 as you have your code set up like that.

请记住，在梯形法则中，您将有限间隔细分为n个部分.第i个 定义为:

Remember, in the trapezoidal rule, you are subdividing the finite interval into n pieces. The i^th piece is defined as:

x_i = a + dx*i;  ,

其中i从1到N-1.请注意，这是从 1 开始的，而不是3.之所以如此，是因为f(a)已经考虑了第一部分，并且我们最多只能将N-1计为N，因为由f(b)占.对于等式，这从2到N，并且通过以此方式修改代码，这正是我们最终要做的.

where i goes from 1 up to N-1. Note that this starts at 1 and not 3. The reason why is because the first piece is already taken into account by f(a), and we only count up to N-1 as piece N is accounted by f(b). For the equation, this goes from 2 to N and by modifying the code this way, this is precisely what we are doing in the end.

因此，您的陈述实际上需要为:

Therefore, your statement actually needs to be:

trapezoidsum = (dx/2)*(f(a) + (2*sum(f(a+dx*(1:n-1))))+f(b));

尝试一下，让我知道您是否获得正确答案. FWIW，MATLAB已经通过使用 trapz 作为@来实现梯形积分ADonda已经指出.但是，在进行此设置之前，您需要正确构造x和y的值.换句话说，您需要先设置dx，然后使用我在上面指定的x_i公式计算x点，然后使用这些点生成y值.然后，使用trapz计算面积.换句话说:

Try this and let me know if you get the right answer. FWIW, MATLAB already implements trapezoidal integration by doing trapz as @ADonda already pointed out. However, you need to properly structure what your x and y values are before you set this up. In other words, you would need to set up your dx before hand, then calculate your x points using the x_i equation that I specified above, then use these to generate your y values. You then use trapz to calculate the area. In other words:

dx = (b-a) / n;
x = a + dx*(0:n);
y = f(x);
trapezoidsum = trapz(x,y);

您可以使用上面的代码作为参考，以查看您是否正确地实现了梯形规则.您的实现和使用上面的代码应生成相同的结果.您要做的就是更改n的值，然后运行此代码以生成曲线下方不同细分区域的面积近似值.

You can use the above code as a reference to see if you are implementing the trapezoidal rule correctly. Your implementation and using the above code should generate the same results. All you have to do is change the value of n, then run this code to generate the approximation of the area for different subdivisions underneath your curve.

我弄清楚了为什么您的代码无法正常工作.原因如下:

I figured out why your code isn't working. Here are the reasons why:

1. for循环是不必要的.看一下for循环迭代.您有一个从i = [3:n]开始的循环，但您完全没有引用循环中的i变量.因此，您根本不需要这个.

1. The for loop is unnecessary. Take a look at the for loop iteration. You have a loop going from i = [3:n] yet you don't reference the i variable at all in your loop. As such, you don't need this at all.

您没有正确计算连续间隔.您需要做的是计算第n个子间隔的梯形和，然后递增此值n，然后再次计算梯形规则.在您的while循环中，此值未正确增加，这就是您的区域从未改善的原因.

You are not computing successive intervals properly. What you need to do is when you compute the trapezoidal sum for the nth subinterval, you then increment this value of n, then compute the trapezoidal rule again. This value is not being incremented properly in your while loop, which is why your area is never improving.

您需要将前一个区域保存在while循环内，然后在计算下一个区域时，即确定两个区域之间的差异是否小于公差.我们还可以在开始尝试并计算n = 2的区域时删除该代码.不需要，因为我们可以将其放在您的while循环中.这样，您的代码应如下所示:

You need to save the previous area inside the while loop, then when you compute the next area, that's when you determine whether or not the difference between the areas is less than the tolerance. We can also get rid of that code at the beginning that tries and compute the area for n = 2. That's not needed, as we can place this inside your while loop. As such, this is what your code should look like:

---

function t=trapintegral(f,a,b,tol)
format long; %// Got rid of format compact. Useless
%// n starts at 2 - Also removed syms x - Useless statement
n = 2;
newtrap = ((b-a)/2)*(f(a) + f(b)); %// Initialize
oldtrap = 0; %// Initialize to 0
while (abs(newtrap-oldtrap)>=tol)
    oldtrap = newtrap; %//Save the old area from the previous iteration
    dx = (b-a)/n; %//Compute width
    %//Determine sum
    trapezoidsum = (dx/2)*(f(a) + (2*sum(f(a+dx*(1:n-1))))+f(b));
    newtrap = trapezoidsum; % //This is the new sum
    n = n + 1; % //Go to the next value of n
end
t = newtrap;
end

通过运行代码，这就是我得到的:

By running your code, this is what I get:

trapezoidsum = trapintegral(@(x) (x+x.^2).^(1/3),1,4,0.00001)

trapezoidsum = 

6.111776299189033

注意

看看我定义函数的方式.您必须使用逐个元素的操作，因为循环内的sum命令将被矢量化.专门查看^操作.您需要在操作前加一个点.完成此操作后，我将获得正确的答案.

Caveat

Look at the way I defined your function. You must use element-by-element operations as the sum command inside the loop will be vectorized. Take a look at the ^ operations specifically. You need to prepend a dot to the operations. Once you do this, I get the right answer.

您说过要至少一个for循环.这是非常低效的，并且在代码中指定一个for循环的人实际上都不知道MATLAB是如何工作的.但是，您可以使用for循环来累积sum项.因此:

You said you want at least one for loop. This is highly inefficient, and whoever specified having one for loop in the code really doesn't know how MATLAB works. Nevertheless, you can use the for loop to accumulate the sum term. As such:

function t=trapintegral(f,a,b,tol)
format long; %// Got rid of format compact. Useless
%// n starts at 3 - Also removed syms x - Useless statement
n = 3;
%// Compute for n = 2 first, then proceed if we don't get a better
%// difference tolerance
newtrap = ((b-a)/2)*(f(a) + f(b)); %// Initialize
oldtrap = 0; %// Initialize to 0
while (abs(newtrap-oldtrap)>=tol)
    oldtrap = newtrap; %//Save the old area from the previous iteration
    dx = (b-a)/n; %//Compute width
    %//Determine sum
    %// Initialize
    trapezoidsum = (dx/2)*(f(a) + f(b));

    %// Accumulate sum terms
    %// Note that we multiply each term by (dx/2), but because of the
    %// factor of 2 for each of these terms, these cancel and we thus have dx
    for n2 = 1 : n-1
        trapezoidsum = trapezoidsum + dx*f(a + dx*n2);
    end
    newtrap = trapezoidsum; % //This is the new sum
    n = n + 1; % //Go to the next value of n
end
t = newtrap;
end

---

祝你好运！

---

Good luck!

这篇关于通过MATLAB中的梯形和进行积分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！