Численные методы и обработка данных

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

Пример: решить систему линейных уравнений

Решение возможно одним из способов (s₁, s₂, s₃ или s₄ – см. приведенную ниже программу).

--> a=[2, 1, 0, 1; 1, -3, 2, 4; -5, 0, -1, -7; 1, -6, 2, 6]; b=[8 9 -5 0];

--> s1=b/a', s2=a\b', s3=b*a'^(-1), s4=b*inv(a')

s1 =

8.1481481 - 1.5185185 11.703704 - 6.7777778

s2 =

8.1481481

- 1.5185185

11.703704

- 6.7777778

s3 =

8.1481481 - 1.5185185 11.703704 - 6.7777778

s4 =

8.1481481 - 1.5185185 11.703704 - 6.7777778

Решить систему линейных уравнений вида можно с помощью функции linsolve : linsolve(a,b). В нашем примере:

a=[2, 1, 0, 1; 1, -3, 2, 4; -5, 0, -1, -7; 1, -6, 2, 6]; b=[8; 9 ;-5 ;0]; x=linsolve(a,-b)

x =

8.1481481

- 1.5185185

11.703704

- 6.7777778

Если система уравнений имеет бесчисленное множество решений, то выводится одно из них:

a=[2, 1; 4, 2]; b=[-7; -14 ];x=linsolve(a,b)

x =

2.8

1.4

Если система не имеет решений:

a=[2, 1; 4, 2]; b=[-7; -13 ];x=linsolve(a,b)

WARNING:Conflicting linear constraints!

x =

[]

Вычисление корней полинома

Функция roots(c) возвращает вектор-столбец из корней полинома с.

Пример: решить уравнение

--> x=[7, 0, 12, 23]; d=roots(x)

d =

0.5564046 + 1.6257442i

0.5564046 - 1.6257442i

- 1.1128093

Примечание: Коэффициенты полинома следует вводить в порядке убывания степеней переменной x. Если в уравнении отсутствует слагаемое, содержащее, например, x², то в векторе коэффициентов на соответствующем месте надо ввести 0.

Решение нелинейных уравнений вида f(x)=0

Уравнения бывают алгебраическими и трансцендентными. Алгебраическим называют уравнение вида . Если уравнение нельзя свести к алгебраическому заменой переменных, то его называют трансцендентным. Пример:

Для решения уравнений, в том числе трансцендентных, в Scilab применяют функцию fsolve(x0,f)

где x0 - начальное приближение, f - функция, описывающая левую часть уравнения f(x)=0.

Пример: решить уравнение

Набираем в окне редактора файл:

function y=f(x)

y=7*x.^3+45*x.^2+12*x+23;

endfunction

и сохраняем его под именем f.sci.Загружаем его в Scilab(Execute/Load into Scilab).

Для нахождения отрезка [a, b], на котором отделен корень данного уравнения, построим график функции .

-->x=-8:0.1:-5; plot(x, f(x)); xgrid()

Из графика видно, что корень отделен на отрезке [-6.5 , -6]. Найдем его, используя функциюfsolve:

-->x0=-6.5;x1= fsolve(x0,f)

Получаем:

x1 =

- 6.2381997

Систему нелинейных уравнений также можно решить, используя функцию fsolve.

clc

function [y]=ff(x)

y(1)=x(1)^2+x(2)^2-1;

y(2)=x(1)^3-x(2);

endfunction

t=fsolve([-.5,-.5],ff)

t =

- 0.8260314 - 0.5636242

Поиск минимума функции y=f(x) на интервале [a, b]

Функция [f1,xopt]=optim(costf,x0) возвращает локальный минимум функции costf.

Функция возвращает минимум функции (значение f1) и точку, в которой этот минимум достигается (xopt).

Главной особенностью функции optim является структура функции costf, которая должна быть такой:

function [f,g,ind]=costf(x,ind)

f=gg(x);//функция, минимум которой мы ищем

g=numdiff(gg,x); //градиент функции f

Endfunction

Если возвращаемое сформированной функцией costf значение ind равно 2, 3 или 4, то функция costf обеспечивает поиск минимума, т.е. в качестве результата функции optim возвращается f и xopt. Если ind=1, то в функции optim ничего не считается, условие ind<0 означает, что минимум f(x) не может быть оценен, а ind=0 прерывает оптимизацию. Вообще говоря, значение параметра ind является внутренним параметром для связи между optim и costf, для использования optim необходимо помнить, что параметр ind должен быть определен в функции costf.

Пример: найти минимум функции .

Решение. Построим график функции для определения интервалов [a, b], на которых находятся экстремумы этой функции.

-->x=-3:.1:3; y=x.^4-3*x.^2-5*x-4; plot(x, y); xgrid()

Из графика видно, что это отрезок [1, 2]. Набираем в окне редактора и отправляем на выполнение файл

function [f,g,r]=z(x,r)

f=x.^4-3*x.^2-5*x-4

g=4*x.^3-6*x-5

endfunction

x0=1;

[fmin,xmin]=optim(z,x0)

Получаем

--> xmin =

1.5233402

fmin =

- 13.193373

Возможен другой вариант, без ручного вычисления производной:

function y=gg(x)

y=x.^4-3*x.^2-5*x-4;

endfunction

function [f,g,r]=z(x,r)

f=gg(x)

g=numdiff(gg,x)

endfunction

x0=1;

[fmin,xmin]=optim(z,x0)

xmin =

1.5233402

fmin =

- 13.193373

В случае функции двух переменных:

(Поиск минимума функции Розенброка )

clc

x0=[-2;2];

function y=gg(x)

y=100*(x(2)-x(1)^2)^2+(1-x(1))^2;

endfunction

function [f,g,r]=z(x,r)

f=gg(x)

g=numdiff(gg,x)

endfunction

[fmin,xmin]=optim(z,x0)

xmin =

0.9999955

0.9999910

fmin =

2.010D-11