Skip to content

Численное дифференцирование в 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 переменных ( в нашем случае переменных больше )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
    
  • Метод нахождения производной четвертого порядка точности для второй производной

  • Формула:4-ый порядок вторая производная

    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».

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

Опубликовано вMatlab

Будьте первым, кто оставит комментарий

    Добавить комментарий

    Ваш e-mail не будет опубликован. Обязательные поля помечены *