Численное дифференцирование
% Работа стандартной функции diff при использовании точных данных
h =.001; x = 0:h:pi; y=sin(x.^2); dy=diff(y)/h;
% Приближенная производная dy.
n=length(x); X=x(1:n-1); DY=2*cos(X.^2).*X; % Истинная производная
figure(1); plot (X,DY, X(1:20:end),dy(1:20:end), '.')
title('Точная и разностная производная')
%
% Возмутим значения функции случайной ошибкой с уровнем delta:
%
delta=0.001
Y=sin(x.^2)+delta*randn(size(x)); dy=diff(Y)/h;
%
% Посмотрим на все это:
%
figure(2); subplot(2,1,1); plot(x,y, X(1:20:end), Y(1:20:end),'.')
title ('Точная и приближенная функция')
legend('Точная','Поиближ.',3)
subplot (2,1,2); plot(X,dy,'.-', X,DY,'y')
title('Точная и приближенная производная')
legend('Приближ.','Точная',3)
nev=norm(DY-dy)/norm(DY)
pause
%
% Попробуем продифференцировать сплайн, построенный по
% возмущенным данным на разных сетках
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
close all;
kk=0; Err=zeros(1,10);
for k=[2 10 30 70 100]
H=k*0.001; kk=kk+1;
xx=0:H:pi; yy=sin(xx.^2)+delta*randn(size(xx));
YY=interp1(xx,yy,x,'spline'); dY=diff(YY)/h;
Err(kk)=norm(dY-DY)/norm(DY);% Ошибка дифференцирования
% Err(kk)=norm(dY(1:end-100)-DY(1:end-100))/norm(DY(1:end-100));
figure; plot(X,dy,'r.-',X,dY,'y',X,DY,'k')
title('Точная и приближенные производные')
legend('Разностная','Сплайновая','Точная',3);
text(1,6,'Шаг ='); text(1.4,6,num2str(H));
pause(1)
end
% Найдем ошибку дифференцирования
figure; plot(0.001*[2 10 30 70 100], Err(1:kk),'o-');
title ('Зависимость ошибки дифференцирования от шага');
xlabel('Шаг'); ylabel('Ошибка дифференцирования')
pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Есть ли наилучший шаг?
% Можно ли и как использовать сплайн для устойчивого
% дифференцирования? Ваше мнение? Попробуйте проделать
% дифференцирование не сплайна, а линейного интерполянта,
% полученного с помощью interpl ('linear'). Где получилось лучше и
% почему?
% Найдите зависимость ошибки дифференцирования от шага в этом случае.
% ОБРАТИТЕ ВНИМАНИЕ: здесь нужно вычислять ошибку по-другому
% Err(kk)=norm(dY{1:end-100)-DY(1:end-100))/norm{DY{lrend-100));
% Это связано с работой interpl. Строка такого вида уже запасена.
% Раскомментируйте ее в выделенном выше фрагменте!
%
close all
%Вычисление второй производной
%
h =.001; x = 0:h:pi; y=sin(x.^2); n=length(x);
d2y=(y(3:n)-2*y(2:n-1)+y(1:n-2))/h^2; % Приближенная производная
X=x(1:n-2); D2Y=2*cos(X.^2)-4*(X.^2).*sin(X.^2); % Истинная производная
figure(1); plot (X,D2Y, X(1:20:end), d2y(1:20:end), '.')
title('Точная и разностная вторая производная')
% Возмутим данные дифференцирования
delta=0.00001
Y=sin(x.^2)+delta*randn(size(x));
d2y=(Y(3:n)-2*Y(2:n-1)+Y(1:n-2))/h^2;
figure(2); subplot(2,1,1); plot(x,y, x(1:20:end),Y(1:20:end),'.')
title('Точная и приближенная функция')
legend('Точная','Приближ.',3)
subplot(2,1,2); plot(X(1:N-2),d2y,'.-', X,D2Y,'y')
title('Точная и приближенная вторая производная')
legend('Приближ.','Точная',2)
nev2=norm(D2Y-d2y)/norm(D2Y)
% Вычисление какой производной более устойчиво: первой или второй?
Приближенные методы решения систем дифференциальных уравнений
Примеры численного решения дифференциальных уравнений
Во многих случаях бывает нужно исследовать так называемые динамические системы - системы, движение которых определяется дифференциальными уравнениями. Рассмотрим простой пример — движение частицы в поле сил, причем ограничимся движением вдоль одной прямой. Координата частицы x и ее скорость v определяются уравнениями
где a(x) = F(x)/m, F(x) — действующая на частицу сила, m — масса частицы. Задача состоит в вычислении зависимостей x (t), v (t) при условии, что заданы начальные значения координаты и скорости х (0) и v( 0).
Вид этих уравнений представляет собой очевидный намек на возможный способ вычислений: выбрать достаточно малое конечное значение величины dt, а затем воспользоваться соотношениями
Многократно применяя эти соотношения, можно рассчитать значения х и v в ряде дискретных, но достаточно близких друг к другу точек. Чтобы оценить величину отклонения от истинного закона движения, запишем, например для х, более точное
равенство
Из него видно, что неточность на каждом шаге пропорциональна (dt)2; для продвижения на время t необходимо число шагов, равное t/dt, так что неточность окажется пропорциональна t × dt. Если необходимо (при заданном интервале t) улучшить точность в 10 раз, то нужно в 10 раз уменьшить шаг dt. Во столько же раз вырастет число шагов и время, необходимое для расчета.
Смысл равенства (*) состоит в учете ускорения на интервале dt. Можно добиться примерно той же точности, что и в (*), если взять среднее значение скорости на том же интервале или, что удобнее всего, значение скорости в середине интервала, v { t + dt/2). Ускорение же, необходимое для выполнения сдвига на шаг по времени для скорости, нужно будет вычислять в середине интервала (t + dt/2, t + dt/2 + dt), т.е. в момент t + dt, что нас тоже устраивает, так как это позволит найти x(t + dt).
Таким образом, для увеличения точности достаточно вычислять значения координат в моменты времени
t t + dt, t + 2dt, t + 3dt,...,
а значения скорости — в моменты
t + dt/2, t + 3dt/2, t + 5 dt/2, t + 7dt/2,...,
Для оценки неточности расчета запишем (вновь только для х)
что совпадает с (*). Неточность имеет порядок (dt)3.
Очевидно, что такой способ расчетов применим и для векторных величин.
Однако этот прием не срабатывает, если сила, действующая на частицу, зависит не только от координат, но и от скорости.