Лабы по Matlab 1- 10 / Zanyatie_7_Reshenie_tipovykh_zadach_algebry_i_anal
.pdfЗанятие 7. Решение типовых задач алгебры и анализа
MATLAB обладает большим набором специальных функций, реализующих различные численные методы. Нахождение корней уравнений и систем, приближение табличных функций, интегрирование, решение задач линейной алгебры, оптимизация, решение систем обыкновенных дифференциальных уравнений – вот далеко не полный перечень возможностей, предоставляемых MATLAB. Мы рассмотрим решение типовых задач на применение вычислительных методов.
Задачи линейной алгебры
В системе MATLAB для решения систем линейных уравнений предусмотрены знаки операций / и \. Чтобы решить систему линейных уравнений:
a11x1 +... +a1n xn =b1
...
an1x1 +... + ann xn =bn
которую обычно кратко записывают в виде:
A·X = B
где A – заданная квадратная матрица размером NxN,a B – заданный вектор – столбец длины N, достаточно применить операцию \ и вычислить выражение A\B.
Операция \ называется левым делением матриц и, будучи примененная к матрицам A и B в виде A\B, эквивалентна вычислению выражения inv(A)*B. Здесь под inv(A) понимается вычисление матрицы, обратной к матрице A.
Для вычисления определителя квадратной матрицы предназначена встроенная функция det(A).
Собственные числаλi и собственные векторыui ≠ 0 квадратной матрицы А удовлетворяют равенствам A ui = λiui . Функция eig, вызванная с входным аргументом матрицей, находит все собственные числа матрицы и записывает их в выходной аргумент – вектор:
>>A = [2 3; 3 5];
>>lam = eig(A) lam =
0.1459
6.8541
Для одновременного вычисления всех собственных векторов и чисел следует вызвать eig с двумя выходными аргументами.
>> [U, Lam] = eig(A) U =
-0.8507 0.5257 0.5257 0.8507
Lam =
0.1459 0
0 6.8541
Первый выходной аргумент и представляет собой матрицу, столбцы которой являются собственными векторами. Для доступа, например, к первому собственному вектору следует использовать индексацию при помощи двоеточия – U(:,1).
Вторым выходным аргументом Lam возвращается диагональная матрица, содержащая собственные числа исходной матрицы.
Проверьте, правильно ли найдены, например, второе собственное число и соответствующий ему собственный вектор. Воспользуйтесь определением:
>> A*U(:, 2) - Lam(2, 2)*U(:, 2) ans =
1.0e-015 *
0.4441 -0.8882
Нахождение нулей функций
Решение уравнения F(x)=0, или нахождение нулей функции, осуществляется с помощью функции fzero(name,x0). В качестве первого аргумента ей передается имя функции, задающей исходное уравнение, вторым аргументом служит начальное приближение к корню. Возвращаемым значением функции fzero является нуль функции name в окрестности точки x0. Определим, в качестве примера, нули функции cos(x) на отрезке от 0 до pi. В качестве начального приближения примем x0=1.
>> x=fzero('cos', 1) x =
1.5708
Если требуется найти корень функции, отличной от стандартной (встроенной в систему MATLAB) и тем самым не имеющей в рамках системы MATLAB фиксированного имени, то можно создать собственную файл-функцию.
Пусть, например, требуется найти корни уравнения cos(x)=x , что эквивалентно нахождению нулей функции, вычисляемой по формуле y=cos(x)−x, не имеющей в рамках системы MATLAB фиксированного имени. В этом случае, прежде всего, нужно создать в редакторе M-файлов функцию, набрав две строки следующего кода:
function y=MyFunction1(x) y=cos(x)-x;
и запомнить их в файле MyFunction1.m. После этого можно воспользоваться функцией fzero:
>> x=fzero('MyFunction1',pi/2) x =
0.7391
Обратите внимание, что имя функции name, задающее исходное выражения, для fzero (это также справедливо для большинства других функций MATLAB, таких как fplot, fminbnd и др.) можно задавать как:
•как имя m-файла (в апострофах и без расширения m) – 'MyFunction1', 'cos'
•как указатель на функцию – @MyFunction1, @cos
•или как формула с одним неизвестным x (y,a и т.п.) – 'cos(x)-x', 'cos(x)'
Еще раз подчеркнем, что функция fzero находит нули только вещественно значных функций одной вещественной переменной. Однако часто бывает необходимо найти комплексные корни вещественнозначных функций, особенно в случае многочленов. Для этой цели в системе MATLAB существует специальная функция roots, которой в качестве аргумента передается массив коэффициентов многочлена.
Например, для многочлена х4–Зх3+Зх2-Зх+2, имеющего два вещественных (1 и 2) и два комплексных корня (i и -i), нужно сначала сформировать массив его коэффициентов (расположив по порядку убывания степени х):
>> Coef=[1,-3,3,-3,2] |
|
|
||
Coef = |
-3 |
3 |
-3 |
2 |
1 |
||||
после чего вызвать функцию roots:
>> r =roots(Coef) r =
2.0000
-0.0000 + 1.0000i -0.0000 - 1.0000i 1.0000
Функция polyval предназначена для вычисления значения полинома от некоторого аргумента:
>> polyval(Coef,1.25) ans =
-0.4805
Аргумент может быть матрицей или вектором, в этом случае производится поэлементное вычисление значений полинома и результат представляет матрицу или вектор того же размера, что и аргумент. Для проверки правильности работы roots, можно вычислить значение полинома от вектора его корней: polyval(Coef,r)
В задаче о нахождении нулей функции сложным моментом является нахождение начального приближения к нулю функции, а также априорная оценка их количества. Поэтому важно параллельно с применением функций типа roots или fzero визуализировать поведение искомых функций на том или ином отрезке значений аргумента (построить график функции).
Поиск минимума функции
В системе MATLAB имеются специальные функции для поиска минимума заданных функций. При этом возможен поиск минимума как для функции одной вещественной переменной, так и для функций многих переменных.
Для функций одной переменной их минимумы разыскивает функция fminbnd: fminbnd(name,x0,xl)
Здесь name представляет имя функции, у которой находятся минимумы, а x0 и xl задают отрезок поиска. Рассмотрим использование fminbnd на примере функции у=3x4-4x3- 12x2. Для этого сначала построим ее график на отрезке [-2,3]:
>> fplot('3*x^4-4*x^3-12*x^2',[-2,3])
Из рисунка 1 видно, что функция y(x) имеет два минимума на отрезках [-2,0] и [0,3], найдем их значения:
>>x1=fminbnd('3*x^4-4*x^3-12*x^2',-2,0) x1 =
-1.0000
>>x2=fminbnd('3*x^4-4*x^3-12*x^2',0,3) x2 =
2.0000
Рис. 1 График функции у=3x4-4x3-12x2
Для поиска минимума функции нескольких переменных применяется функция fminsearch:
xmin=fminsearch(name,x0)
Здесь name является именем функции нескольких переменных, для которой ищется минимум, а x0 - это вектор ее аргументов, с которого начинается поиск. Для иллюстративного примера создадим простую функцию двух переменных имеющую минимумом точку (0,0).
function y = MyFunc2(x) y = x(1)^2+ x(2)^2;
Этот текст надо записать в файл MyFunc2.m. После этого можно вызвать функцию fminsearch:
>> xmin = fminsearch (@MyFunc2, [1,1]);;
которая приближенно находит вектор xmin координат точки минимума:
>>xmin(1) ans = -2.1024e-005
>>xmin(2) ans =
2.5485e-005
Обе найденные координаты близки к своим точным значениям, равным нулю. Для функций нескольких переменных еще важнее, чем для ранее рассмотренных функций одной вещественной переменной, постараться априорно оценить количество и приблизительное нахождение локальных минимумов, и тут могут существенно помочь трехмерные графики.
Управление ходом вычислений
Функции fzero, fminbnd и fminsearch допускают определение дополнительных параметров для управления вычислительным процессом и контроля за ним. Формат вызова:
fzero(name, x0, options) fminbnd(name, x1, x2, options) fminsearch(name, x0, options)
Параметры задаются в управляющей структуре, которую мы будем называть options, как в справочной системе MATLAB, хотя имя может быть произвольным. Перед вызовом вычислительных функций следует предварительно сформировать переменную options в соответствии с характером требуемого контроля, воспользовавшись функцией optimset. Переменная options на самом деле является структурой. Это новый тип данных – до сих пор вы имели дело только с числовыми массивами. Это аналог типа sruct в C++, и record в Pascal. В рамках данного курса умение работать со структурами не требуется, можно просто следовать приведенным правилам заполнения, просмотра и использования управляющей структуры options.
Приступим к формированию структуры options на примере минимизации функции одной переменной при помощи fminbnd. Для изменения точности следует сформировать структуру options, указав требуемую точность (например 10-9) по аргументу при помощи параметра TolX:
>>format long
>>options = optimset('TolX', 1.0e-09);
>>x2 = fminbnd('exp(x.^2)+sin(3*pi*x)', -0.5, 0, options)
x2 = -0.16289942841268
В общем случае входные аргументы optimset задаются попарно:
options = optimset('Свойство1', Знач1, 'Свойство2', Знач2, ...);
Некоторые возможные сочетания свойств и их значений приведены в следующей таблице:
Свойство |
Значение |
|
Результат |
|
Примечание |
|
Display |
'off' |
Информация о вычислительном |
|
|||
|
|
процессе не выводится |
|
|
||
|
'iter' |
Выводится |
информация |
о |
|
|
|
|
каждом |
шаге |
вычислительного |
|
|
|
|
процесса |
|
|
|
|
|
'final' |
Выводится информация только о |
|
|||
|
|
завершении |
вычислительного |
|
||
|
|
процесса |
|
|
|
|
|
'notify' |
Выводится |
предупреждение, |
|
||
|
|
если |
процесс не сходится |
|
||
|
|
(используется по умолчанию) |
|
|
||
MaxFunEvals |
Целое число |
Максимальное |
количество |
Только для |
||
|
|
вызовов исследуемой функции |
|
fminbnd |
||
|
|
|
|
|
|
fminsearch |
MaxIter |
Целое число |
Максимальное |
количество |
Только для |
|
|
итераций |
вычислительного |
fminbnd |
|
|
процесса |
|
fminsearch |
TolFun |
Положительное |
Точность по |
функции для |
Только для |
|
вещественное |
останова вычислений |
fminbnd |
|
|
число |
|
|
fminsearch |
TolX |
Положительное |
Точность по аргументу функции |
|
|
|
вещественное |
для останова вычислений (по |
|
|
|
число |
умолчанию 10-6) |
|
|
Ограничивать количество вызовов функций и число итераций имеет смысл, если есть опасение, что получить решение не удастся из-за расхождения вычислительного процесса. В некоторых случаях приходится уменьшать точность, например, если вычисление исследуемой функции занимает много времени, а требуется получить только несколько первых значащих цифр ответа.
Вычисление определенных интегралов
Для вычисления интегралов методом трапеций в системе MATLAB предусмотрена функция trapz:
Integ = trapz(х,у);
Одномерный массив х (вектор) содержит дискретные значения аргументов подынтегральной функции. Значения подынтегральной функции в этих точках сосредоточены в одномерном массиве у. Чаще всего для интегрирования выбирают равномерную сетку, то есть значения элементов массива х отстоят друг от друга на одну и ту же величин (шаг интегрирования). Точность вычисления интеграла зависит от величины шага интегрирования: чем меньше этот шаг, тем больше точность. Вычислим простой интеграл методом трапеций с шагом интегрирования рi/10:
∫π cos(x)dx
0
>>dx = pi/10;
>>x = 0:dx:pi;
>>y=cos(x);
>>I1 = trapz(x,y)
I1 =
3.8858e-016
Обычно для достижения высокой точности требуется выполнять интегрирование с очень малыми шагами, а контроль достигнутой точности осуществлять путем сравнения последовательных результатов. При одном и том же шаге интегрирования методы более высоких порядков точности достигают более точных результатов.
В системе MATLAB метод интегрирования более высоких порядков точности реализуются функцией: quad (метод Симпcона). Этот метод являются к тому же адаптивным. Последнее означает, что пользователю нет необходимости контролировать достигнутую точность результата путем сравнения последовательных значений, соответствующих разным шагам интегрирования. Все это функция выполняет самостоятельно.
Как и многие другие функции системы MATLAB, функция quad может принимать различное количество параметров. Минимальный формат вызова этой функции
quad(name, x1, x2)
включает в себя три параметра: имя подынтегральной функции – name, нижний предел интегрирования – x1 и верхний предел интегрирования – x2. Кстати, если эта функция не может обеспечить получение необходимой точности (расходящийся или близкий к этому интеграл), то она возвращает символическую бесконечность Inf.
Для примера, найдем значение интеграла
∫1 e−x sin(x)dx
−1
Первым шагом является создание функции, вычисляющей подынтегральное выражение, ее текст:
function f = fint(x) f = exp(-x).*sin(x);
Теперь для вычисления интеграла вызовите quad, задав первым аргументом ссылку на функцию fint, а вторым и третьим — нижний и верхний пределы интегрирования. В качестве выходного аргумента можно указать имя переменной, в которую следует записать найденное значение:
>>format long
>>I = quad(@fint,-1,1)
I = -0.66349146785310
По умолчанию функция quad вычисляет приближенное значение интеграла с точностью 10-6. Для изменения точности вычислений следует задать дополнительный четвертый аргумент:
>> I = quad(@fint,-1,1,1.0e-07)
I=
-0.66349366574399
Задания для самостоятельной работы
Решить систему уравнений AX=B. Произвести проверку полученного решения. Вычислить определитель, собственные значения и собственные вектора матрицы A. Провести проверку для полученных собственных векторов и собственных значений.
=−9
1.−4x1 − x2 −5x3 = −1x1 +5x2 + x3 = −94x2 −5x33x1 +−
x1 +5x2 − x3 = 0
3.−5x3 = −10
−2x1 + x2 − x3 = 5
Варианты
|
x + 4x |
2 |
− x = −14 |
||
2. |
|
1 |
|
3 |
|
4x1 +5x2 −3x3 = −24 |
|||||
|
x + 2x |
2 |
− x = −8 |
||
|
|
1 |
|
3 |
|
|
4x −5x −5x =5 |
||||
4. |
|
1 |
|
2 |
3 |
−3x1 −3x2 + 4x3 =10 |
|||||
|
−5x + 2x |
−5x = −9 |
|||
|
|
1 |
|
2 |
3 |
|
−5x |
−5x |
+5x = −5 |
|
−3x |
+ x + 2x = −4 |
|||||||
5. |
|
|
1 |
2 |
|
|
3 |
6. |
|
|
1 |
2 |
3 |
−4x1 + x2 = −1 |
− x1 − x2 + 4x3 = −4 |
||||||||||||
|
5x1 + x2 −5x3 =9 |
|
−3x1 −3x2 − x3 = −12 |
||||||||||
|
2x |
+ x + x |
|
= 4 |
|
−2x |
−5x |
+5x =1 |
|||||
7. |
|
1 |
|
2 |
3 |
|
8. |
|
|
1 |
2 |
3 |
|
− x1 +3x2 −3x3 = −3 |
5x1 |
− x2 +5x3 = 7 |
|||||||||||
|
−4x |
+3x |
+3x = −18 |
|
−3x |
+ 2x |
−5x = −5 |
||||||
|
|
|
1 |
2 |
|
|
3 |
|
|
|
1 |
2 |
3 |
|
−5x |
−3x |
−2x = 22 |
|
−3x +5x + x = −8 |
||||||||
9. |
|
|
1 |
2 |
|
|
3 |
10. |
|
|
1 |
2 |
3 |
5x1 |
− x2 +5x3 = −24 |
−2x1 + 4x2 − x3 = −11 |
|||||||||||
|
−3x |
+ 4x |
|
= −5 |
|
4x |
−3x −3x = 2 |
||||||
|
|
|
2 |
3 |
|
|
|
|
1 |
|
2 |
3 |
|
Задания для самостоятельной работы
Найти корни многочлена. Выполнить проверку полученного решения.
Варианты
1. |
p = x8 + 2x7 −3x5 +11x2 −5 |
2. |
p = −5x8 −6x6 − x5 +6x3 + 2x2 −12 |
3. |
p =16x8 − x7 −3x5 +8x + 2 |
4. |
p = −3x8 −3x6 − x4 + x3 + 2x −6 |
5. |
p = 2x7 −3x5 + 4x3 +8x −8 |
6. |
p = 2x7 +6x4 −2x2 +8 |
7. |
p = x9 +9x8 + x3 +12x −1 |
8. |
p = −2x10 −3x6 −4x4 + x3 −2x2 +6 |
9. |
p = 2x8 +15x6 +7x3 −9x −4 |
10. |
p =16x5 −7x4 −9x2 +3 |
Задания для самостоятельной работы
Исследовать поведение функций графически на заданном отрезке. Вычислить определенный интеграл заданных подынтегральных функций по промежутку между:
•соседними минимумами или
•соседними нулями или
•между соседними минимумом и нулем
(в зависимости от поведения функции) методом трапеций и методом Симсона. Все вычисления вести с точностью 10-10.
Варианты
1. |
y = sin x cos x |
x [0;2π] |
|
2. |
y = sin2 x cos2 x |
−0.1 x [0;2π] |
||||||||
|
x2 +1 |
|
|
|
|
|
x3 +1 |
|
|
|
||||
3. |
y = ln(x +1) |
ex +e−x |
x [−3;3] |
|
4. |
y = −x4 − x3 +6x2 +0.5 |
x [−1;2.2] |
|||||||
5. |
y = sin x − x2 cos x x [−5;−0.5] |
|
6. |
y = −1 −sin2 x + x2 cos x |
x [−5;5] |
|||||||||
7. |
y = ex 2 +sin3πx x [0;1.5] |
|
8. |
y = ex 2 −sin |
πx |
+10x3 x [−2.2;2] |
||||||||
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
9. |
y = 2 − |
|
1 |
|
x [0;3] |
10. |
y = |
1 |
|
+ |
|
1 |
|
−3 x [0;2.5] |
(x −0.3)2 +0.1 |
|
(x −0.3)2 +0.1 |
|
(x −0.9)2 +0.04 |
||||||||||