Лекции по Matlab

function y=f1(z)

  y=z*exp(-z)+sin(z) 

     Для нахождения не только корня, но и значения функции в нем, можно обратиться к fzero с 2 выходными параметрами:

     >> [x, f]=fzero('x.*exp(-x)+sin(x)',[3, 4])

     x =

         3.2665

     f =

       2.0817e-016 

     Иногда, например, при нахождении корня уравнения tg(x)=0 функция fzero выдает неожиданный результат:

     >> [x, f]=fzero('tan(x)',[1, 2])

     x =

         1.5708

     f =

      -1.2093e+015

     Судя  по значению функции найденный x не является корнем. 

     Полный  формат использования функции fzero:

     [x, f, e_flag, inform]=fzero(fun, x0)

Если  e_flag = 1, то удалось найти интервал, на концах которого функция fun меняет знак и e_flag=-1 в противном случае. Структура inform содержит 3 поля с информацией о числе итераций, выполненных при поиске корня, количестве обращений к функции fun, алгоритме, использованном для решения.

     >> [x, f, e_flag, inform]=fzero('x.*exp(-x)+sin(x)',[3, 4])

     x =

         3.2665

     f =

       2.0817e-016

     e_flag =

          1

     inform =

         iterations: 8

          funcCount: 8

          algorithm: 'bisection, interpolation' 

     Для алгоритма, примененного в fzero принципиальным является принятие функций значений разных знаков на концах интервала. Например, уравнение x²=0 не может быть решено с помощью fzero (если не указать в качестве приближения явно 0). 

     Решение систем нелинейных уравнений 

     Для решения систем нелинейных уравнений  вида F(x)=0 предназначена функция fsolve. Здесь x – вектор или матрица неизвестных, F – функция, значением которой является вектор или матрица.

     Алгоритм  работы fsolve использует начальное приближение x0 и базируется на минимизации суммы квадратов компонент функции F методами Гаусса-Ньютона и Левенберга-Марквардта.

     Простейшее  использование: x=fsolve(F,x0). В частности fsolve можно использовать вместо функции fzero для нахождения корня единственного нелинейного уравнения:

     >> x=fsolve('sin(x)',0.1)

     В MATLAB существует 2 версии данной функции, поэтому при подобном кратком вызове fsolve будет выведна на несколько строчек пояснительная информция с предложением выбора одной из версий функции с помощью более точного указания аргументов (для поиска корня в этом случае будет использована более свежая версия fsolve). Для этого необходимо использовать 3 параметр:

     x=fsolve('sin(x)',0.1,optimset('Display','off'))

     x =

       1.2495e-011 

     С помощью функции optimset задаются свойства (параметры), влияющие на ход итерационных процессов (оптимизации, минимизации). Задание значения любого из свойств производится парой параметров, первый – наименование, второй – значение:

     options=optimset('Name1',Val1,'Name2',Val2,…);

     x=fxxx(fun, x0,…, options);

     В вышеописанном примере был установлен режим отображения информации о  шагах итерационного процесса (Display). Он может быть off (вывод отключен), iter (вывод на каждой итерации), notify (вывод в случае необходимости), final (финальный вывод). 

     Рассмотрим  пример решения следующей системы  уравнений:

     Объявим в m-файле функцию funcsc, значения которой сформируем в виде вектора-столбца:

     function y=funcsc(x)

       y=[x(1)+x(2)-sin(pi*x(1));x(1)-x(2)-cos(pi*x(2))];

     Далее вызов fsolve

     >> [x, f]=fsolve(@funcsc,[0.2;1],optimset('Display','off'))

     x =

         0.3915

         0.5510

     f =

       1.0e-005 *

         0.9147

        -0.0339

     Если  вызвать fsolve с нулевым приближением -1 и -2, то, судя по возвращаемому значению f, процедура поиска попала в локальный минимум:

     >> [x, f]=fsolve(@funcsc,[-1;-2],optimset('Display','notify'))

     Optimizer is stuck at a minimum that is not a root

     Try again with a new starting guess

     x =

        -0.4595

        -1.6778

     f =

        -1.1454

         0.6883 

 

     Численное дифференцирование  функции одной  переменной 

     Производная – понятие локальное. Поэтому, если функция задана таблицей значений, для нахождения ее производной в каком-либо из узлов обычно используют интерполяционные многочлены невысоких степеней, хорошо приближающие функцию в окрестности этого узла. Для построения таких многочленов используют формулы, основанные на таблицах разностей. Наиболее известны интерполяционные формулы Ньютона.

     Конечные  разности

     Пусть на оси абсцисс отмечены узлы x₁, x₂, …, x_n и функция задана таблицей значений y_i=f(x_i).

     Таблицы разностей:

     Если  узлы равноотстоящие, то составляются также разности 2, 3 и т.д. порядков:

 

     В MATLAB функция diff строит по заданному вектору x вектор разностей указанного порядка:

     dx=diff(x) %разности 1-го порядка

     dx=diff(x,n) %разности n-го порядка

     Например, для функции y=log(x) составим разности первого и второго порядка:

     >> x=1:0.1:2;

     >> dx=diff(x);

     >> y=log(x);

     >> dy=diff(y);

     >> d2y=diff(y,2); 

     Пусть узлы в которых заданы значения функции  равноотстоящие: Δx₁=Δx₂=…=Δx_n=h. Тогда можно построить интерполяционные формулы, позволяющие приближенно находить значения функции вблизи произвольного узла x_i (q=( x- x_i)/h).

     Формулы Ньютона 1 порядка (линейное приближение):

     y=f(x)~y_i+q*Δy_i (интерполяция вперед)

     y=f(x)~y_i+q*Δy_i-1 (интерполяция назад)

     Формула Ньютона 2 порядка (квадратичное приближение):

     y=f(x)~y_i+q*Δy_i+(q*(q-1)*Δ²y_i)/2! (интерполяция вперед)

     y=f(x)~y_i+q*Δy_i-1+(q*(q+1)*Δ²y_i-2)/2! (интерполяция назад)

     Формул  Стирлинга 1 порядка (линейное приближение):

     y=f(x)~y_i+q*(Δy_i-1 + Δy_i)/2

     Формул  Стирлинга 2 порядка (квадратичное приближение):

     y=f(x)~y_i+q*(Δy_i-1 + Δy_i)/2 + q²*Δ²y_i-1/2!

     Если  продифференцировать полусумму формул Ньютона 1-го порядка для интерполяции вперед и назад или формулу Стирлинга 1-го порядка, то получим:

     f'(x_i)= (Δy_i-1 + Δy_i)/(2h)= (y_i+1 - y_i-1)/(2h)

     (Для  концевых узлов i=1 и i=n данная формула неприменима, в этом случае приходится применять по отдельности формулы Ньютона, желательно более высокого порядка:

     f'(x_i)= (Δy₁ - Δ²y₁/2)/h= (-3y₁ + 4y₂ - y₃)/(2h)

     f'(x_n)= (Δy_n-1 - Δ²y_n-2/2)/h= (3y_n - 4y_n-1 + y_n-2)/(2h) 

     Задача: Вычислить значения производной функции log(x) в узлах сетки аналитически и по вышеприведенным формулам.

 

     Численное интегрирование 

     В основе всех способо численного интегрирования  (вычисления площади криволинейной  трапеции, расположенной под графиком заданной функции) лежит суммирование значений функции в некоторых узлах, причем значения берутся с определенными весами, зависящими от положения узлов на отрезке интегрирования.

     Простейшая  формула численного интегрирования с равноотстоящими узлами – формула  трапеций:

     

     Она является точной для линейных функций (из вида остаточного члена).

     Для вычисления интеграла по этой формуле  в MATLAB используется функция trapz.

     >> x=1:0.1:2;

     >> y=log(x);

     >> trapz(x,y)

     ans =

         0.3859

     Для сравнения точное значение интеграла  равно 0.3863.

     Для функции trapz узлы могут быть и не равноотстоящими: