В системе MATLAB имеется функция diff, предназначенная для вычисления конечных разностей; ее синтаксис:
z=diff(y), где – массив; функция diff возвращает конечные разности смежных элементов y:
;
число элементов на единицу меньше числа элементов х. Для перехода к производной конечные разности необходимо разделить на шаг независимой переменной х;
z=diff(y,n) – вычисление конечных разностей n-го порядка
Следует иметь в виду, что процедуре численного дифференцирования присуща методическая погрешность, возрастающая при увеличении шага независимой переменной и порядка производной. Поэтому функцию diff следует применять лишь для приближенного вычисления производных невысокого порядка относительно медленно изменяющихся (гладких) зависимостей.
Составим программу (deriv.m), иллюстрирующую вычисление производных до третьего порядка включительно:
n=input('n ==> '); % Ввод высшего порядка производной
x=-3:step:(3+n*step); % Учет уменьшения числа элементов
% при дифференцировании
m=length(x)-n; % Число элементов массива высшей
% производной
y=exp(-x.^2); % Дифференцируемая функция
z1=diff(y)/step; % Первая производная
if n==2
z2=diff(z1)/step; % Вторая производная
End
if n==3
z2=diff(z1)/step;
z3=diff(z2)/step; % Третья производная
z2=z2(1:m);
End
x=x(1:m); % Уменьшение размерности массивов
y=y(1:m);
z1=z1(1:m);
figure % Открываем графическое окно
if n==1 % Строим графики исходной функции
plot(x,y,x,z1) % и ее производных
End
if n==2, plot(x,y,x,z1,x,z2), end
if n==3, plot(x,y,x,z1,x,z2,x,z3,'k'), end
Grid
Title('Численное дифференцирование')
Xlabel('Число элементов массивов')
Ylabel('Функция и ее производные')
Pause
Close
Упрощенный вариант программы (deriv1.m) выводит графики как зависимости от номеров элементов. В этом случае корректировать их число с учетом уменьшения при дифференцировании не нужно.
step=input('step ==> ');
n=input('n ==> ');
x=-3:step:3;
y=exp(-x.^2);
color=['r','g','k'];
Figure
Plot(y), hold on
for k=1:n
z=diff(y,k)/step.^k;
Plot(z,color(k)), hold on
End
Grid
Title('Численное дифференцирование')
legend('y(x)','dy/dx','d^2y/dx^2','d^3y/dx^3')
Xlabel('Число элементов массивов')
Ylabel('Функция и ее производные')
Pause
Close
Упражнение. Составить программу дифференцирования другой функции и проанализировать влияние величины шага переменной на точность вычисления производной.
Дифференцирование может быть организовано и в командном режиме. Если задано аналитическое выражение исходной функции, во многих случаях операцию дифференцирования можно реализовать точно – в символьной форме.
5.4.2. Функции eval, feval
Возможность интерпретации и вычисления выражений, заданных в символьной форме, придает системе MATLAB дополнительную гибкость.
Функция eval(‘Выражение’)возвращает результат выполнения выражения, заданного символьно, при известном значении независимой переменной
script %eva
s='exp(-0.5*x).*sin(3*x)';
x=0:0.1:10;
y=eval(s);
Plot(x,y)
Grid
Pause
Close
Функция feval(‘Выражение’,x)возвращает значение внешней функции,заданной в символьной форме, при аргументе (массиве) х. Эту функцию применяют, когда внешняя функция вычисляется в теле другой процедуры, например, при интегрировании дифференциальных уравнений. В следующем простом примере использована созданная ранее внешняя функция func.
script %feva
Global nn
nn=2;
f=’func’;
x=0:0.1:4;
Plot(x,feval(f,x))
Grid
Pause
Close
5.4.3. Решение обыкновенных дифференциальных уравнений (ОДУ)
В инженерной практике необходимость решения ОДУ возникает часто. Система MATLAB имеет встроенный решатель, который позволяет выбирать конкретный способ решения, задавать начальные условия (задача Коши) и оптимизировать вычислительный процесс. Решатель также используется в пакете SimuLink. Кроме того, при инсталляции системы MATLAB 5.х можно подключить специализированный пакет для решения дифференциальных уравнений (ДУ) в частных производных.
Вычислительный алгоритм при интегрировании ОДУ любого порядка заключается в том, что это уравнение превращают в систему ДУ первого порядка с соответствующими правыми частями: , где y –вектор переменных состояния; t – независимая переменная (в инженерных задачах чаще всего время); – вектор правых частей.
Решатель ОДУ использует различные численные методы, названия (имена) которых должны заменять слово solver:
где результат вычислений [t,y] – матрица, столбцы которой представляют собой массивы t, а также и т.д.; ‘f’ – имя внешней функции, вычисляющей правые части ДУ первого порядка; задает диапазон изменения независимой переменной (от t0до tk); y0 – вектор начальных значений; options – аргумент, создаваемый с помощью функции odeset.
В MATLAB реализуются следующие методы численного интегрирования:
Гладкие системы:
ode45– алгоритм Рунге-Кутта 4-го и 5-го порядков; рекомендуется для пробного интегрирования системы ДУ;
ode23– алгоритм Рунге-Кутта 2-го и 3-го порядков; применяется при низких требованиях к точности;
ode113– решатель переменного порядка, основанный на формуле Адамса; эффективнее ode45 при высоких требованиях к точности и сложных выражениях в правых частях системы ДУ;
Жесткие системы: ode15s, ode23s, ode23t, ode23tb.
В качестве примера приведем программу для решения ДУ Ван-дер-Поля . . Введем обозначения: ; тогда исходное уравнение можно представить как систему из ДУ первого порядка:
Сначала составляем внешнюю программу (van.m) для вычисления правых частей ДУ:
% Вычисление правых частей
function ypr=van(t,y)
ypr=[y(2,:); (1-y(1,:).^2)*y(2,:)-y(1,:)];
затем script-файл (vander.m) (начальные условия: при ):
% Интегрирование ДУ Ван-дер-Поля
script %vander
y0=[0 1]';
[t,y]=ode45('van',[0 30],y0);
Figure
Plot(t,y(:,1),t,y(:,2))
title(‘Решение уравнения Ван-дер-Поля’)
xlabel(‘t’), ylabel(‘y, dy/dt’)
legend(‘ y’,’ dy/dt’)
Grid
Pause
Plot(y(:,1),y(:,2))
Grid
title(‘Фазовый портрет колебаний’)
xlabel(‘y’), ylabel(‘dy/dt’)
Close
Рекомендуется изучать тексты встроенных m-файлов (они выводятся по команде type name; кроме того, можно вызывать тексты и напрямую, загружая их в окно редактора). Это – образцы тщательно разработанных фирмой MathWorks программ с использованием всех возможностей системы MATLAB.
СИМВОЛЬНЫЕ ВЫЧИСЛЕНИЯ
Символьные, т.е. аналитические, вычисления в среде MATLAB 5.x реализуются, если инсталлирован пакет Symbolic Math Toolbox(его файлы размещаются в папке toolbox\symbolic). В этом случае используется до 100 наиболее часто применяемых функций символьной системы Maple V R4, и MATLAB становится универсальной системой, возможности которой расширяются еще более благодаря прямому доступу к ядру Maple. Ниже рассматриваются возможности версии пакета Symbolic 2.0.1.
Для получения справки по пакету необходимо набрать команду » help symbolic.Перед началом работы рекомендуется вызвать следующие демонстрационные программы: symintro – начальное ознакомление с пакетом Symbolic; symcalcdemo–исчисление; symlindemo–символьная линейная алгебра; symeqndemo– символьное решение уравнений.
Рассмотрим наиболее важные варианты символьных вычислений в системе MATLAB.