Составное правило Симпсона

У меня есть этот код для составного правила Симпсона. Тем не менее, я довольно долго возился с ним, и, похоже, я не могу заставить его работать.

Как исправить этот алгоритм?

function out = Sc2(func,a,b,N)
% Sc(func,a,b,N)
% This function calculates the integral of func on the interval [a,b]
% using the Composite Simpson's rule with N subintervals.
x=linspace(a,b,N+1);
% Partition [a,b] into N subintervals
fx=func(x);
h=(b-a)/(2*N);
%define for odd and even sums
sum_even = 0;
for i = 1:N-1
   x(i) = a + (2*i-2)*h;
   sum_even = sum_even + func(x(i));
end

sum_odd = 0;

for i = 1:N+1
   x(i) = a + (2*i-1)*h;
   sum_odd = sum_odd + func(x(i));
end
% Define the length of a subinterval
out=(h/3)*(fx(1)+ 2*sum_even + 4*sum_odd +fx(end));
% Apply the composite Simpsons rule
end

matlab numerical-methods numerical-analysis numerical-integration
person Jake    schedule 12.11.2014    source источник
comment
Кажется, я не могу заставить его работать, и как я могу исправить этот алгоритм? не являются конкретными или полезными на этом или другом сайте? Покажите, как вы тестируете эту функцию, т.е. какую функцию вы интегрируете и как вы вызываете свою Sc2 функцию. Если вы получаете сообщения об ошибках, сообщите о них полностью. Как вы определяете, верен ли ваш результат? Отредактируйте свой вопрос. И если это домашнее задание, укажите это.   -  person horchler    schedule 13.11.2014

Ответы (1)


arrow_upward
6
arrow_downward

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

Кроме того, вы суммируете от 1 до N+1 или от 1 до N-1 для нечетных или четных значений. Это неверно. Помните, что вы выбираете любое другое значение в нечетном или четном интервале, поэтому на самом деле это должно происходить в цикле от 1 до N/2 - 1. Чтобы не выяснять, на что умножать i, просто пропустите это и сделайте свой цикл пошаговым по 2.

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

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

Таким образом, вот модифицированная версия вашей функции с оптимизацией, которую я имею в виду:

function out = Sc2(func, a, b, N)

h = (b – a) / N; %// Width of each interval
odd = 1 : 2 : n-1; %// Define odd interval 
xodd = a + h*odd; %// Create odd x values 
even = 2 : 2 : n-2; %// Create even interval
xeven = a + h*even; % Create even x values 

%// Return area
out = (h/3)*(func(a) + 4*sum(func(xodd)) + 2*sum(func(xeven))+ func(b));

Однако, если вы хотите, чтобы ваш код работал, вам просто нужно изменить пределы итераций for цикла, а также значение h. Вам также необходимо удалить некоторые строки кода и изменить имена некоторых переменных. Следовательно:

function out = Sc2(func,a,b,N)
% Sc(func,a,b,N)
% This function calculates the integral of func on the interval [a,b]
% using the Composite Simpson's rule with N subintervals.

%// Define width of each segment
h = (b - a) / N; %// Change

%//define for odd and even sums
sum_even = 0;
for i = 2 : 2 : N-2 %// Change 
   x = a + i*h; %// Change
   sum_even = sum_even + func(x);
end

sum_odd = 0;

for i = 1 : 2 : N-1 %// Change
   x = a + i*h %// Change
   sum_odd = sum_odd + func(x);
end

%// Output area
out = (h / 3)*(func(a) + 2*sum_even + 4*sum_odd + func(b)); %// Change

end
person rayryeng    schedule 13.11.2014