3358

2)Зададим функции f(x,y) и g(x,y).

>f:=(x,y)->evalf(x^7-5*x^2*y^4+1510);

>g:=(x,y)->evalf(y^5-3*x^4*y-105);

3)Для отделения корней системы построим графики неявных функций f(x,y)=0 и g(x,y)=0.

>implicitplot({f(x,y)=0,g(x,y)=0},x=-5..5,y=-5..5);

Рис. 1.4

Видим пять пересечений графиков функций, значит, система имеет, как минимум, пять решений.

4)Определим якобиан

>J:=[[diff(f(x,y),x),diff(f(x,y),y)],

[diff(g(x,y),x),diff(g(x,y),y)]];

5)Введём функцию Jn(x,y), являющуюся определителем матрицы J.

>Jn:=unapply(det(J),x,y);

6)Определим величины An и Bn в виде соответствующих определителей, при этом для исключения символьных преобразований результат оформляем через evalf

>An:=det([[f(x,y),diff(f(x,y),y)],

[g(x,y),diff(g(x,y),y)]]);

>Bn:=det([[diff(f(x,y),x),f(x,y)],

[diff(g(x,y),x),g(x,y)]]);

7)Присвоим начальные значения переменным. x0 и y0 задают точку начального приближения (определяется графически).

20

>x0:=-3.; y0:=-1.; eps:=0.001: n:=0:

8)Организуем цикл вычислений по методу Ньютона

>while abs(evalf(f(x0,y0)))+abs(evalf(g(x0,y0)))>eps do

if abs(evalf(Jn(x0,y0)))<10^(-4) then

print(`Якобиан = 0`); break fi; # Если якоби-

ан = 0, выйти из программы x0:=evalf(x0-subs(x=x0,y=y0,An)/Jn(x0,y0)); y0:=evalf(y0-subs(x=x0,y=y0,Bn)/Jn(x0,y0)); if n=50 then break fi; # Защита от зацикливания n:=n+1;

od:

9) Вывод результата: (x0,y0) – решение системы, n – число итераций

> evalf({x0,y0});n;

{–2.844483289, –.5348543094} 4

10) Для вычисления остальных четырех корней системы нужно выполнить шаги 7–9, предварительно изменив в п. 7 начальное приближение. В частности, установив

> x0:=-2.; y0:=-3.;

получим второе решение системы: {–2.573256593, –2.304767678}

4

Третье, четвертое и пятое решение получим с начальны-

ми приближениями а) x0:=2., y0:=–3., б) x0:=-2., y0:=3.,

в) x0:=2., y0:=3. соответственно. Эти значения следует последовательно подставить в п. 7 и выполнить все последующие команды программы.

В завершение сведем все найденные решения в ответе.

Ответ: (–2,844; –0,535), (–2,573, –2,305), (2,430; –2,873), (–1,922, 2,957), (1,994; 3,012).

Из других методов решения систем нелинейных уравнений имеют применение метод простой итерации, метод наискорейшего спуска, градиентные методы.

21

1.5. Использование стандартных функций системы Maple

Для получения численного решения нелинейного уравнения или системы уравнений удобно использовать Maple-

функцию fsolve(eqns,vars,options);

Здесь eqn, eqn1, … – уравнения, содержащие неизвестные переменные var, var1, … Эта функция может быть использована со следующими параметрами:

complex – находит один или все корни (чаще полинома) в комплексной форме;

fulldigits – определяет счет для полного числа цифр, заданного функцией Digits;

maxsols=n – задает вычисление только n корней;

interval – задается в виде a..b или x=a..b, или {x=a..b, y=c..d,…} и обеспечивает поиск корней в указанном интервале.

Следующие примеры демонстрируют возможности использования fsolve.

> fsolve(exp(x)=sin(x),x=–4..0); # Поиск корня урав-

нения ex sinx, принадлежащего отрезку [–4, 0] –3.183063012

> fsolve(exp(x)=sin(x),x=–7..–4); # То же на отрез-

ке [–7, –4]

–6.281314366

Обратите внимание, Maple может выдать только один корень, лежащий на заданном отрезке, даже если на самом деле их там несколько. Поэтому вся ответственность по отделению корней возлагается на пользователя.

>fsolve(exp(x)=sin(x),x,complex); # Комплексный корень

.3627020561 – 1.133745919 I

>fsolve({sin(x+y)-exp(x)*y=0,x^2-y=2},{x,y},

{x=-1..1,y=-2..0}); # Численное решение системы с указанием интервалов, которым принадлежат неизвестные

{x = –.6687012050, y = –1.552838698}

22

Упражнения

1.Решить нелинейные уравнения и системы уравнений одним из итерационных методов или с помощью встроенной Maple-функции. Предварительно отделить корни графическим методом.

а) Уравнения с одним неизвестным

2. Найти корни нелинейных уравнений, приведенных в таблице, методами

а) половинного деления, б) простой итерации, в) касательных, г) секущих.

Количество и положение корней определить графически. Оценить точность полученныхзначений.

23

24

3. Решить систему двух нелинейных уравнений методом Ньютона. Начальные приближения определить графически.

25

4*. Составить программу решения уравнения f(x)=0 с помощью метода:

5*. Для функции f(x) найти точку минимума x* и минимальное значение f(x*) методом Ньютона

26

2. РЕШЕНИЕ ЗАДАЧ ЛИНЕЙНОЙ АЛГЕБРЫ

Методам линейной алгебры принадлежит особая роль в математике, особенно в численном анализе. Ведь многие задачи математического анализа, дифференциальных и интегральных уравнений после дискретизации сводятся к решению задач линейной алгебры. Более того, большинство нелинейных моделей в малой окрестности решения могут быть описаны как конечномерные линейные модели, в основе исследования которых также лежит линейная алгебра.

2.1. Матричные и векторные операции

Прежде чем приступить к решению задач в пакете Maple, напомним некоторые определения, относящиеся к матрицам и векторам.

Матрица размера (m n) – прямоугольная двумерная таблица, содержащая m строк и n столбцов элементов. В расширенной трактовке каждый элемент матрицы может быть не только числом, но и переменной, функцией, выражением и т.д. Квадратная матрица – матрица, у которой число строк равно числу столбцов.

Вырожденная или особенная матрица – квадратная матрица с определителем, равным нулю.

Единичная матрица E – квадратная матрица, у которой диагональные элементы равны 1, а остальные равны нулю.

Нижняя (верхняя) треугольная матрица – матрица, у ко-

торой элементы над (под) главной диагональю равны нулю. Транспонированная матрица AT – матрица A, у которой

столбцы и строки поменялись местами, т.е. элементы матриц AT и A связаны соотношением aijT aji .

Обратная к A матрица – это матрица A–1, которая, будучи умноженной на квадратную матрицу A, дает единичную матрицу E.

След матрицы – сумма её элементов на главной диагонали. Ранг матрицы – наибольший из порядков отличных от

нуля миноров квадратной матрицы.

27

Определитель матрицы – это многочлен от элементов квадратной матрицы A порядка n, определяемый формулой

n

detA a1 j ( 1)j 1M1 j ,

j 1

где M1 j – определитель матрицы порядка n–1, полученной из матрицы A вычеркиванием первой строки и j-го столбца.

Симметричная матрица – матрица A, отвечающая условию AT=A.

Положительно определенная матрица – квадратная мат-

рица, отвечающая условию xAxT >0 для любого вектора x 0. Собственный вектор (n n)-матрицы A – любой вектор x 0 размерности n, удовлетворяющий условию Ax= x, где

– некоторое число, называемое собственным значением (числом) матрицы A. Собственные значения матрицы являются корнями её характеристического многочлена det(A– E).

Норма – обобщенное понятие абсолютной величины числа. В качестве нормы вектора x = (x1, x2, …, xn) чаще всего используют следующие

странстве совпадает с длиной вектора);

n

б) x1 | xi |;

i 1

в) x maxi | xi |.

Нормой матрицы A, обозначаемой ||A||, называется наименьшее из чисел c в неравенстве ||Ax|| c||x|| x 0. В зависимости от применяемой нормы вектора будут получаться различные значения нормы матрицы:

где (ATA) – максимальное собственное число матрицы ATA.

28

Число обусловленности матрицы A определяется какAA 1 и характеризует степень зависимости погрешно-

сти решения системы Ax=b от погрешности правой части. Если матрица близка к особенной, то число будет для такой матрицы велико. В этом случае говорят, что матрица A плохо обусловлена.

>restart; # Очистка оперативной памяти

Вядре системы Maple содержатся весьма небольшие и минимально необходимые средства для решения задач линейной алгебры. Полный набор таких средств реализован в подключаемых пакетах. Основным из них является пакет расширения linalg, который содержит свыше ста функций. Чтобы его подключить, следует запустить команду:

>with(linalg);

Ниже приводятся примеры выполнения основных действий над матрицами и векторами.

1.Задание квадратной матрицы 3-го порядка:

>A:=matrix([[2,1,3],[5,1,0],[7,8,9]]);

2.Вычисление определителя матрицы:

>det(A);

72

3. Формирование матрицы, полученной вычеркиванием i-й строки и j-го столбца матрицы A:

> minor(A,2,1);

Определитель этой матрицы дает минор элемента с индексами i и j матрицы A.

4.Вычисление ранга матрицы:

>rank(A);

3

29