Примеры численного решения дифференциальных уравнений

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

 

% Работа стандартной функции 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)²; для продвижения на время 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)³.

Очевидно, что такой способ расчетов применим и для векторных величин.

Однако этот прием не срабатывает, если сила, действующая на частицу, зависит не только от координат, но и от скорости.