Численное дифференцирование в Matlab

Доброго времени суток! Сегодня мы продолжаем говорить об обработке статических данных (ранее разбирали апроксимацию и интерполяцию в Matlab).

Численное дифференцирование в Matlab. Теория

Процесс дифференцирования подразумевает нахождение производной. Для получения разных порядков производной, можно провести операцию дифференцирования несколько раз. Более детально, о смысле и значение производной, предлагаю ознакомиться здесь. Численное дифференцирование строится на использовании аппарата конечных разностей и соответствующего многообразия аппроксимаций.

Нахождение первой производной в Matlab

Сперва скажу, что есть стандартные средства дифференцирования в Matlab и именно методы (своего рода алгоритмы). К стандартным относится такой инструмент как diff(x). Она принимает аргумент x, который является массивом. Подробнее о diff() читайте здесь.

Методы численного дифференцирования применяются, если исходную функцию y(x) трудно или невозможно продифференцировать аналитически. Методы имеют разную погрешность при расчётах. Рассмотрим основные из них и оценим погрешность помощью Matlab:
Сравним результаты, дифференцируя функцию sin(x).

h=0.2; %определим шаг сетки
x=0:h:pi; %интервал значений х
n=length(x); %число необходимых итераций
dy=cos(x); %производная y = sin(x) => y' = cos(x), так будем судить об отклонении и погрешности.
dz=sin(x); % для сравнения

• Метод нахождения производной правой конечной разностью

Формула:

for i=1:(n-1)
 dy1(i)=(sin(x(i+1))-sin(x(i)))/h;
 er1(i)=abs(dy(i)-dy1(i));
end

• Метод нахождения производной левой конечной разностью

Формула:

for i=2:n
 dy2(i)=(sin(x(i))-sin(x(i-1)))/h;
 er2(i)=abs(dy(i)-dy2(i));
end

• Метод нахождения производной центральной конечной разностью

Формула:

for i=2:(n-1)
 dy3(i)=(sin(x(i+1))-sin(x(i-1)))/(2*h);
 er3(i)=abs(dy(i)-dy3(i));
end

• Метод нахождения производной четвертого порядка точности

Формула для 4-го порядка точности и 4 переменных ( в нашем случае переменных больше )

for i=3:(n-1)
 dy4(i)=(-sin(x(i+1))+27*sin(x(i))-...
 27*sin(x(i-1))+sin(x(i-2)))/(24*h);
 er4(i)=abs(cos(x(i)-0.5*h)-dy4(i)); % абсолютную погрешность, которая вычисляется в точке x(i)-0.5*h
end

%рисуем все на графиках
 plot(x([1:(n-1)]),er1([1:(n-1)]),'-o',...
             x([2:n]),er2([2:n]),'-p',...
     x([2:(n-1)]),er3([2:(n-1)]),'-h',...
     x([3:(n-1)]),er4([3:(n-1)]),'-*');
 title('Погрешность ("разность" анлитического и численного решения)');
 legend('Правая', ' Левая', 'Центральная', '4-ый порядок')
 grid on;
figure; plot(x([1:(n-1)]),dy([1:(n-1)]),'-d',...
                x([1:(n-1)]),dz([1:(n-1)]),'-<',...
                x([1:(n-1)]),dy1([1:(n-1)]),'-o',...
     x([3:(n-1)]),dy4([3:(n-1)]) - 0.5*h,'-*');
 legend('cos(x)', ' sin(x)', 'производная по правой конечной разности', 'производная по 4-ому порядоку')
 grid on

Наши выводы:

Нахождение второй производной в Matlab

Теперь провернём похожие фокусы с использование формулы для центральной конечной разности и формулы 4-го порядка точности для второй производной.
Сравним результаты, дифференцируя функцию sin(x).

h=0.2;  %определим шаг сетки
x=0:h:pi;  %определим сетку
n=length(x);  %число итераций
d2y=-sin(x);   %точное значение второй производной d2y=-sin(x)

• Метод нахождения производной центральной конечной разностью для второй производной

Формула:

for i=2:(n-1)
    d2y1(i)=(sin(x(i+1))-2*sin(x(i))+sin(x(i-1)))/h^2;
    er1(i)=abs(d2y(i)-d2y1(i));<18pan>
end

• Метод нахождения производной четвертого порядка точности для второй производной

Формула:

for i=3:(n-2)
    d2y2(i)=(-sin(x(i+2))+16*sin(x(i+1))-30*sin(x(i))+...
                    16*sin(x(i-1))-sin(x(i-2)))/(12*h^2);
    er2(i)=abs(d2y(i)-d2y2(i));
end

Наши выводы:

Метод Рунге: уточнения формул численного дифференцирования

Я выделил в отдельности данный метод, так как он служит только для уточнения результатов. Быстренько пробежимся по примеру программы:

%Теория и практика в среде MATLAB
%Программа иллюстрирующая метод Рунге по повышению точности численного
%дифференцирования путем комбинирования пары расчетов на двух равномерных сетках,
%причем вторая содержит вдвое большее количество узлов
d=sin(x);
h=0.1;  %определяем шаг исходной сетки
x=0:h:pi;  %формируем исходную сетку
n=length(x);  %количество итераций
hm=h/2; %определяем шаг более подробной сетки
xm=0:hm:pi;  %создаем сетку вдвое более подробную
m=length(xm); %количество итераций для подробной
%рассчитываем производную с помощью правой разности на исходной сетке и оцениваем
for i=1:(n-1)
    dy(i)=(sin(x(i+1))-sin(x(i)))/h;
    er1(i)=abs(cos(x(i))-dy(i));  %соответствующую абсолютную ошибку 
end
%рассчитываем производную с помощью правой разности на более частой сетке и оцениваем
for i=1:(m-1)
    dym(i)=(sin(xm(i+1))-sin(xm(i)))/hm;
    er2(i)=abs(cos(xm(i))-dym(i));
end
%уточняем значение производной с помощью метода Рунге
for i=1:(n-1)
    dyrunge(i)=dy(i)-2*(dy(i)-dym(2*i-1));
    er3(i)=abs(cos(x(i))-dyrunge(i));
end
%строим общий график со всеми тремя кривыми ошибок
plot(x([1:(n-1)]),er1([1:(n-1)]),'-o',...
     xm([1:2:(m-1)]),er2([1:2:(m-1)]),'-p',...
     x([1:(n-1)]),er3([1:(n-1)]),'-h');
  title('Погрешность ("разность" анлитического и численного решения)');
 legend( 'Правая конечная на исходной', 'Правая конечная на частой', 'уточнение Рунге')
 grid on;
 figure;plot(x([1:(n-1)]),d([1:(n-1)]),'-d',...
     x([1:(n-1)]),dy([1:(n-1)]),'-o',...
     xm([1:2:(m-1)]),dym([1:2:(m-1)]),'-p',...
     x([1:(n-1)]),dyrunge([1:(n-1)]),'-h');
 legend('sin(x)', 'Правая конечная на исходной', 'Правая конечная на частой', 'уточнение Рунге')
 grid on;

И наш вывод:

Численное дифференцирование в Matlab. Пример задачи

Теперь рассмотрим задачу.
1. Найти 1-ую и 2-ую производные функции, заданной таблично.

xi -3      -2      -1      0       1       2       3
yi      -0,71   -0,01   0,51    0,82    0,88    0,51    0,49

2. Найти 1-ую и 2-ую производные функции y = ln(cos(x)) в точке х= 0,5 различными методами, а именно рассмотреть формулы простые и многоточечные. Оценить точность аппроксимации (это разница между значением производной, вычисленным по точной формуле, полученной аналитически, и её значением, вычисленным по конечно-разностным формулам).
Исследовать влияние величины шага на точность вычисления производных по различным формулам.

Решаем первый пункт.
Нам не дана сама функция, а только точки через которые она проходит. Это значит, что мы должны найти такую функцию, которая приближенна к данным точкам. Найдём её с помощью апроксимации. Тему апроксимации в Matlab мы разбирали ранее. Выглядеть решение будет так:

x=[-3 -2 -1 0 1 2 3];
y=[-0.71  -0.01  0.51  0.82  0.88  0.51 0.49];

nm=4;   %задаем степень аппроксимирующего полинома 
p=polyfit(x,y,nm);   %полиномизируем данные
y1=polyval(p,x);   %%расчитываем значения y1 в зависимости от коэффициентов полинома
dyt = polyder(p);  %берём производную по полиному
dyt2 = polyder(dyt);  %берём вторую производную
dy = polyval(dyt,x);   
dy2 = polyval(dyt2,x);
plot(x,y,'-*'); %строим точки
hold on
plot(x,y1,'-<');  %апроксимация
plot(x,dy,'-or');  %производная первого порядка
plot(x,dy2,'-og');  %производная второго порядка
grid on
legend('данные точки', 'апроксимирующая функция', 'производная первого порядка', 'производная второго порядка');

Наш вывод:

Если вам необходимо вывести значение переменных в процессе решения, достаточно просто убрать «;» в конце строчки кода.

Решаем второй пункт.

h = 0.2;
x = 0:h:pi;
n = length(x);
y = log(cos(x));
dy = -tan(x);  %аналитическое диф уравнение  
for i=1:(n-1)  % производная по правой конечной разности
    dy1(i) = (log(cos(x(i + 1))) - log(cos(x(i))))/h;
    er1(i) = abs(dy(i) - dy1(i));
end
for i=2:n  % производная по левой конечной разности
    dy2(i) = (log(cos(x(i))) - log(cos(x(i - 1))))/h;
    er2(i) = abs(dy(i) - dy2(i));
end
for i = 2:(n-1) % производная по центральной разности
    dy3(i) = (log(cos(x(i+1)))-log(cos(x(i-1))))/(2*h);
    er3(i) = abs(dy(i)-dy3(i));
end
x(i) - 0.5 * h;
for i = 3:(n-1)  % производная по 4ому порядку
    dy4(i) = (-log(cos(x(i+1)))+27*log(cos(x(i)))-27*log(cos(x(i-1)))...
        +log(cos((x(i-2)))))/(24*h);
    er4(i) = abs(-tan(x(i)-0.5*h)-dy4(i));
end
plot(x(1:(n-1)),er1(1:(n-1)),'-o',...
             x(2:n),er2(2:n),'-p',...
     x(2:(n-1)),er3(2:(n-1)),'-h',...
     x(3:(n-1)),er4(3:(n-1)),'-*');
  legend('Правая', ' Левая', 'Центральная', '4-ый порядок')
 %вторая производная 
 grid on;
 figure;
 d2y = 1\(cos(x).^2);
for i = 2:(n-1)  % производная по центральной разности для второй производной
    d2y1(i) = (log(cos(x(i+1)))-2*log(cos(x(i)))+log(cos(x(i-1))))/(h.^2);
    er1(i) = abs(d2y(i)-d2y1(i));
end
for i = 3:(n-2)  % производная по 4ому порядку
    d2y2(i) = (log(cos(x(i+2)))+16*log(cos(x(i+1)))-30*log(cos(x(i)))+...
        16*log(cos(x(i-1)))-log(cos(x(i-2))))/(12*h.^2);
    er2(i) = abs(d2y(i)-d2y2(i));
end
plot(x(2:(n-1)),er1(2:(n-1)),'-o',...
     x(3:(n-2)),er2(3:(n-2)),'-p');
 legend( 'производная по центральной конечной разности', 'производная по 4-ому порядоку')
 grid on;

В теории такие алгоритмы мы уже подробно разобрали. Этот момент «Исследовать влияние величины шага на точность вычисления производных по различным формулам.» я не расписывал в решении, так как его мы разобрали в теории в Методе Рунге.
Наш вывод:

На этом мы закончим. Если остались вопросы, задавайте их в комментариях. Также вы можете скачать исходники чтобы лучше понять тему: «Численное дифференцирование Matlab».

Скачать исходник первого примера
Скачать исходник второго примера
Скачать исходник третьего примера
Скачать исходник первого пункта
Скачать исходник второго пункта