ИЭ (станции+реле лектор Горячевский) / Лабы / лаба 8 / lab_8

%вычисление определённых интегралов с помощью
%метода левых, правых ,средних прямоугольников,
%метода трапеции и метода Симпсона
function [] = lab_8()

clc;clear; %очистить экран

%по аналогии с предыдущей работой
%в качестве исходного данного возьмем модуль напряжений
%с предыдущих расчетов

U = [
9.94;
10.09;
9.96;
10.01;
10.09;
10.09 ];


f = 50; %частота основной гармоники
w = 2 * pi * f; %круговая частота 314 Гц
T = 1 / f; %период частоты 50 Гц = 20 мс
dt = 0.005; %шаг дискретизации по оси времени


%задаём значения по оси времени как простой массив значений
t = 0:dt:T;
%генерируем функцию в каждой точке оси времени
%значения предыдущих расчетов теперь используются просто как
%коэффициенты
u = U(1)+U(2)*cos(w*t)+U(3)*sin(w*t)+U(4)*cos(2*w*t)+U(5)*sin(3*w*t)+U(6)*cos(4*w*t)

%вычисляем интеграл функции по методу левых прямоугольников
%для этого берем первое значение и умножаем его на ширину шага и т.д.
sum = 0;
for i = 1:length(t)-1
sum = sum + u(i) * dt;
end
U_mlp = sum;

%по методу правых прямоугольников
%действия аналогичные, но на первом шаге берём не первое, а второе
%значение и т.д. (то есть значение справа от текущего шага)
sum = 0;
for i = 1:length(t)-1
sum = sum + u(i+1) * dt;
end
U_mpp = sum;

%для метода средних прямоугольников нужно вычислить значение ф-ции
%которое находится посередине шага. Для этого подставляем его в
%исходную формулу

sum = 0;
for i = 1:length(t)-1

tau = ( t(i)+t(i+1) )/2;
v = U(1)+U(2)*cos(w*tau)+U(3)*sin(w*tau)+U(4)*cos(2*w*tau)+U(5)*sin(3*w*tau)+U(6)*cos(4*w*tau);
sum = sum + v * dt;
end
U_msp = sum;

%метод трапеции оказывается проще по записи и точнее, поэтому
%используют в основном его
sum = 0;
for i = 1:length(t)-1
sum = sum + (u(i)+u(i+1))/2 * dt;
end
U_trap = sum;

%формула Симпсона является сочетанием метода трапеций и средних
%прямоугольников
sum = 0;
for i = 1:length(t)-1

tau = ( t(i)+t(i+1) )/2;
v = U(1)+U(2)*cos(w*tau)+U(3)*sin(w*tau)+U(4)*cos(2*w*tau)+U(5)*sin(3*w*tau)+U(6)*cos(4*w*tau);
sum = sum + 1/6 *( u(i)+4*v+u(i+1) ) * dt;
end
U_sim = sum;

disp('МЛП:')
U_mlp
disp('МПП:')
U_mpp
disp('МСП:')
U_msp
disp('МТ:')
U_trap
disp('МС:')
U_sim
end

%Отчёт загрузить до 18:00 10.12.2021
%Попробовать разный шаг дискретизации от 0.02 до 0.0001