Причина, по которой этот код не работает, заключается в том, что в вашем суммировании для правила трапеций есть две небольшие ошибки. Я имею в виду именно это утверждение:

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

Вспомним уравнение трапецеидального правила интегрирования:

^{Источник: Википедия}

Для первой ошибки f(x) должно быть f(a), поскольку вы включаете начальную точку, и его не следует оставлять символическим. На самом деле вам следует просто избавиться от оператора syms x, так как он бесполезен в вашем скрипте. a соответствует x1 в соответствии с приведенным выше уравнением.

Следующая ошибка — второй член. На самом деле вам нужно умножить значения индекса (3:n-1) на dx. Кроме того, на самом деле это должно идти от (1:n-1), и я объясню позже. Приведенное выше уравнение переходит от 2 к N, но для наших целей мы собираемся перейти от 1 к N-1, поскольку ваш код настроен таким образом.

Помните, что в правиле трапеций вы делите конечный интервал на n частей. i^ая часть определяется как:

x_i = a + dx*i;  ,

где i идет от 1 до N-1. Обратите внимание, что это начинается с 1, а не с 3. Причина в том, что первая часть уже учитывается f(a), и мы считаем только до N-1, так как часть N учитывается f(b). Для уравнения это изменяется от 2 до N, и, изменив код таким образом, именно это мы и делаем в конце.

Таким образом, ваше заявление на самом деле должно быть:

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

Попробуйте это и дайте мне знать, если вы получите правильный ответ. FWIW, MATLAB уже реализует трапециевидную интеграцию, выполняя trapz как @ADonda уже указал. Однако вам необходимо правильно структурировать значения x и y, прежде чем настраивать это. Другими словами, вам нужно будет заранее настроить dx, затем рассчитать x очков, используя уравнение x_i, которое я указал выше, а затем использовать их для получения значений y. Затем вы используете trapz для вычисления площади. Другими словами:

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

Вы можете использовать приведенный выше код в качестве справки, чтобы увидеть, правильно ли вы реализуете правило трапеций. Ваша реализация и использование приведенного выше кода должны давать одинаковые результаты. Все, что вам нужно сделать, это изменить значение n, а затем запустить этот код, чтобы сгенерировать аппроксимацию области для различных подразделений под вашей кривой.

Редактировать - 17 августа 2014 г.

Я понял, почему ваш код не работает. Вот почему:

1. Цикл for не нужен. Взгляните на итерацию цикла for. У вас есть цикл, идущий от i = [3:n], но вы не вообще ссылаетесь на переменную i в своем цикле. Соответственно, вам это совершенно не нужно.

2. Вы неправильно вычисляете последовательные интервалы. Что вам нужно сделать, так это, когда вы вычисляете трапецеидальную сумму для nth подинтервала, вы затем увеличиваете это значение n, а затем снова вычисляете трапецеидальное правило. Это значение не увеличивается должным образом в вашем цикле while, поэтому ваша область никогда не улучшается.

3. Вам нужно сохранить предыдущую область внутри цикла while, а затем, когда вы вычисляете следующую площадь, именно тогда вы определяете, меньше ли разница между областями, чем допуск. Мы также можем избавиться от этого кода в начале, который пытается вычислить площадь для n = 2. Это не нужно, так как мы можем поместить это в ваш цикл while. Таким образом, вот как должен выглядеть ваш код:

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

Запустив ваш код, я получаю следующее:

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

trapezoidsum = 

6.111776299189033

Предостережение

Посмотрите, как я определил вашу функцию. Вы должны использовать поэлементные операции, так как команда sum внутри цикла будет векторизована. Обратите особое внимание на операции ^. Вам нужно поставить точку перед операциями. Как только вы сделаете это, я получу правильный ответ.

Редактировать № 2 — 18 августа 2014 г.

Вы сказали, что хотите хотя бы один цикл for. Это очень неэффективно, и тот, кто указал наличие в коде одного цикла for, на самом деле не знает, как работает MATLAB. Тем не менее, вы можете использовать цикл for для накопления термина sum. В качестве таких:

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

Удачи!