Нелинейные уравнения и системы, интегрирование, решение обыкновенных дифференциальных уравнений и систем, обработка результатов эксперимента

Нелинейные уравнения и системы. Численное решение уравнения f(x)=0 часто проводят в два этапа. На первом этапе отделяют корни уравнения, то есть определяют интервалы изоляции корня. Например, график функции пересекает ось Х и окрестности точки пересечения (или нескольких таких точек) это и есть интервалы изоляции корня. На втором этапе проводят уточнение, то есть сужают интервал изоляции корня до ширины, равной заданной погрешности по Х.

Алгебраическое уравнение, оно же полином, оно же многочлен a₀xⁿ+a₁x^n-1+a₂x^n-2+…+a_n-1x+a_n=0

Его вектор коэффициентов p=[a₀ a₁ a₂…a_n-1 a_n]

Произведение двух полиномов вычисляет функция conv(p1,p2) которой передаются векторы коэффициентов обоих полиномов. Деление полиномов (первого на второй) осуществляет функция [q r]=deconv(p1,p2), которая возвращает ответ и остаток от деления.

Если надо вычислить значение полинома при определенном значении х, то для этого используется функция polyval(p1,x). Нахождение производной от полинома осуществляется функцией polyder.

Решить алгебраическое уравнение вида a₀xⁿ+a₁x^n-1+…+a_n-1x+a_n=0 можно функцией x=roots(p); где x вектор корней алгебраического уравнения (полинома), p вектор коэффициентов полинома, то есть р=[a₀ a₁ …a_n-1 a_n]. Например, для уравнения х²-1=0 вектор р=1 0 -1 тогда х=roots(p)  и получается, что х=1 -1. Это и есть корни уравнения.

Например, для нахождения корней многочлена функцией roots(p), пусть вектор p3 = 1 0 0 1. Наоборот, построить полином по вектору его корней может функция poly(x).

function Start_p1p2 p1=[2 0 1]; p2=[1 0 0 1]; p=conv(p1,p2) p3=deconv(p,p1) a1=polyval(p3,2) a2=polyder(p3) x=roots(p3) a3=poly(x) end p = 2 0 1 2 0 1

p3 = 1 0 0 1 a1 = 9 a2 = 3 0 0 x = -1.0000          0.5000 + 0.8660i 0.5000 - 0.8660i a3 = 1.0000 -0.0000 0.0000 1.0000

Уравнение, в котором неизвестное входит в аргумент трансцендентных функций, называется трансцендентным уравнением. Для решения их, предварительно приведенных к виду f(x)=0, используют функцию fzero. Формат ее вызова таков: [x,y]=fzero(name,x0); где name имя М-функции, вычисляющей левую часть уравнения вида f(x)=0.

x0 начальное приближение или интервал его изоляции [a, b];

Например пусть f(x)= 0.2e^x -2(х-1)² тогда см. рис.1.25

На рис.1.25 показано, что при трех разных интервалах изоляции корня найдены три различных корня, которые были помещены в массив х по одному, по мере их нахождения, как х(1), х(2) и х(3), а затем выведены в строчку все три как содержимое массива х.

Аналогично поступили со значениями функции у(1), у(2) и у(3). Их нахождение было важно для контроля того, что действительно значения функции очень близки к нулю (они все получились порядка 10^-15, что практически и есть ноль, с учетом вычислительных погрешностей).

Слева на рис.1.25 видно окно Workspace, где отображены значения переменных и их тип. В принципе, в данной задаче эта информация дублируется в обоих окнах (то есть в Command Window и в Workspace).

Рис. 1.25. Решение уравнения при различных интервалах изоляции корня

Содержимое функции func продемонстрировано на рис.1.26.

Рис. 1.2 6. Функция func(x)

Рекомендуется по возможности построить график. Потом уже описать функцию в М-файле. Затем использовать функцию fzero для нахождения корня.

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

При вызове функции fzero возвращаются только вещественные числа. Комплексные (мнимые) корни она не вычисляет.

function Start_exFzero1 function y=exFzero_1(x) y=((x-1).^2).^(1/3)-(x.^2).^(1/3); end x=[-2:0.1:2]; y=((x-1).^2).^(1/3)-(x.^2).^(1/3); plot(x,y,'r-'); grid on; fzero(‘exFzero_1’,[0 1] end
В Command Window >> ans=     0.5000

Рис.1.27. График трансцендентного уравнения и решение уравнения функцией fzero.

С помощью функции polyval можно получить графическое (рис.1.28) решение уравнения.

function exPolyval p=[1 0.4 0.6 -1.6]; x=[-1:0.1:3]; y=polyval(p,x); plot(x,y,'k-'); grid on; end

Рис. 1.28. Функция polyval обеспечивает графическое решение уравнения.

Численное интегрирование. Для интегрирования существует очень много алгоритмов в современной вычислительной математике. Будем рассматривать выборочно некоторые из них.

Интегрирование методом трапеций означает, что подинтегральная площадка разбивается на трапециевидные подплощадки, границы между трапециями это вертикальные линии, проведенные через x ₁, x ₂, …, x_n равноотстоящие точки на оси Х, идущие от начала до конца интервала интегрирования [ x ₁ x_n ]. Верхние границы трапеций это отрезки, соединяющие точки f (x ₁), f (x ₂), …, f (x_n). Нижние границы трапеций ось Х.

Для метода трапеций есть функция trapz (x, y).

При интегрировании методом Симпсона верхние границы подплощадок являются отрезками параболы. Для него есть функция quad (‘ name ’, a, b, tol, trace), где перечислены аргументы имя М-функции, левая и правая границы интегрирования, точность вычислений, а также параметр задающий таблицу с промежуточными точками результатами вычислений. Есть также весьма точный алгоритм Ньютона-Котеса 8-го порядка. Он вызывается функцией quad 8(‘ name ’, a, b, tol, trace).

function Start_Integr x=[0:0.1:5]; function y=exInt(x) y=x.^2+x-3; end y=exInt(x); a1=quad('exInt',0,5, 0e-06,5) I=trapz(x,y) end   9 0.0000000000 1.35790000e+000 -2.3171466272 11 1.3579000000 2.28420000e+000 14.1273143810 13 3.6421000000 1.35790000e+000 27.3564989128

a1 = 39.1667 I = 39.1750

Пример решения определенного интеграла функцией trapz (методом трапеций) для вышеописанной функции func (см. рис.1.26), где для обработки вектора аргументов поставлена точка перед знаком возведения в степень y=0.2*exp(x) -2*(x-1).^2  при пределах интегрирования [0 1] с шагом разбиения 0,1 приведен на рис.1.29. Вполне очевидно, что с ростом числа точек разбиения по х (числа шагов, то есть при уменьшении шага по х) будет расти точность вычисления интеграла.

Рис. 1.29. Интегрирование методом трапеций.

На рис.1.29 наглядно показано разбиение интервала интегрирования по х с заданным шагом, а также вектор значений функции в точках разбиения по х (то есть при данных значениях аргумента). Получилось значение интеграла равное -0,3261.