Решение систем нелинейных уравнений.
Для решения системы нелинейных уравнений вида F(x)=0, где F(x) – векторная функция, а x – радиус вектор можно использовать функцию fsolve из пакета Optimization Toolbox. Рассмотрим в качестве примера систему из двух уравнений:
,
.
x0 = [-5; -5]; % Нулевое приближение
F ='[2*x(1)-x(2)-exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];'; % Вектор-функция
options=optimset('Display','final'); % Задание управляющих параметров
x = fsolve(F, x0, options) % Вызов оптимизирующей функции
Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun
x =
0.5671
0.5671
Преобразование Фурье.
Дискретное преобразование Фурье (DFT) является основным средством цифровой обработки сигнала. MATLAB содержит функции быстрого преобразования Фурье (FFT) для одномерных, двумерных и многомерных массивов. Прямое и обратное преобразования Фурье основаны на формулах:
,
,
.
Прямое преобразование Фурье (функция fft) переводит описание сигнала (функции времени) из временной области в частотную, а обратное преобразование Фурье (функция ifft) – из частотной области во временную.
Создадим, например, вектор времени и сигнал:
t = (0:1/99:1); % Вектор времени
x = sin(2*pi*15*t) + sin(2*pi*40*t); % Сигнал
Найдем частотное описание сигнала, его амплитуду и фазу:
y = fft(x); % Вычисление DFT of x
m = abs(y); p = unwrap(angle(y)); % Расчет амплитуды и фазы
Построим графики зависимости амплитуды и фазы от частоты:
subplot(1,2,1);
f = (0:length(y)-1)*99/length(y); % Частотный вектор
plot(f, m); title('Magnitude');
set(gca,'XTick',[15 40 60 85]);
subplot(1, 2, 2);
plot(f,p*180/pi); title('Phase');
set(gca,'XTick',[15 40 60 85]);
Рассмотрим применение преобразования Фурье для фильтрации сигнала. Создадим двухчастотный сигнал в 50Гц и 120Гц на фоне сильного шума, создаваемого генератором случайных чисел:
t = 0:0.001:0.6;
x = sin(2*pi*50*t)+sin(2*pi*120*t);
y = x + 2*randn(size(t));
Для выделения исходного сигнала переведем его в частотное представление:
Y = fft(y, 512);
Pyy = Y.* conj(Y) / 512;
Здесь использован n-точечный формат вызова функции fft. Метод является особенно эффективным при n=2m, где m – целое положительное число.
Построим графики зашумленного и восстановленного сигналов:
subplot(1, 2, 1);
plot(y(1:50))
title('Искаженный сигнал')
Xlabel('time (seconds)')
subplot(1, 2, 2);
f = 1000*(0:256)/512;
plot(f, Pyy(1:257))
title('Плотность спектра')
xlabel('частота (Hz)')
На временном графике полезный сигнал полностью забит шумом, тогда как спектрограмма сигнала показывает два явных пика на общем шумовом фоне.
Решение обыкновенных дифференциальных уравнений.
Для решения систем обыкновенных дифференциальных уравнений (ОДУ) в MATLAB реализованы различные численные методы – решатели (solvers) ОДУ. Рассмотрим вначале задачу Коши, записанную в векторной форме:
Перечислим решатели предназначенные для нежестких систем ОДУ:
ode45 – одношаговые явные методы Рунге-Кутта 4-го и 5-го порядка точности. Это классический метод, рекомендуемый для начальной пробы решения;
ode23 – одношаговые явные методы Рунге-Кутта 2-го и 3-го порядка. При низких требованиях к точности этот метод может оказаться более эффективным, чем ode45;
ode113 – многошаговый метод Адамса-Башворта-Мултона переменного порядка. Это адаптивный метод, который может обеспечить высокую точность решения.
Не все трудные задачи определяются жесткими системами, но все жесткие системы являются трудными для неспециализированных решателей. Перечислим решатели жестких систем:
ode15s – многошаговый метод переменного порядка (от 1-го до 5-го) основанный на формулах численного дифференцирования. Это адаптивный метод, его стоит применять, если решатель ode45 дает сбой;
ode15s – одношаговый метод, использующий модифицированную формулу Розенброка 2-го порядка. Обеспечивает высокую скорость вычислений при низкой точности;
ode23t – метод трапеций с интерполяцией. Используется для решения умеренно-жестких задач;
ode23tb – реализация неявного метода Рунге-Кутта. При низкой точности этот метод может оказаться более эффективным, чем ode15s.
Покажем применение решателя ОДУ на классическом примере уравнения Ван-дер-Поля, записанного в виде системы двух дифференциальных уравнений первого порядка:
Опишем правую часть системы уравнений для случая =1 с помощью функции vdp1:
function dydt = vdp1(t, y)
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
Для решения системы на отрезке [0 20] с начальными условиями y(1) = 2 и y(2) = 0 будем использовать решатель ode45:
[t, y] = ode45(@vdp1,[0 20],[2; 0]);
Построим график решения:
plot(t, y(:, 1), '-', t, y(:, 2), '--')
title('Solution of van der Pol Equation, \mu = 1');
xlabel('time t');
ylabel('solution y');
legend('y_1','y_2')
