Как найти корни трансцендентного уравнения в scilab

Основы работы в SciLab. На примере экзаменационных вопросов по КСВЭ

Завтра мне сдавать экзамен по такому предмету как КСВЭ(Компьютерный сервис вычислительного эксперимента). А лучший способ подготовки — это написание статьи. Я рассмотрю часть вопросов к экзамену, которые связаны с SciLab.
Статья больше подходит для студентов, использующих scilab для проверки решения, для или для сдачи экзамена по дисциплине КСВЭ. Для более подробного изучения надо читать литературу, указанную в конце статьи

Основные термины
Scilab (читается Сайлэб) — пакет прикладных математических программ, предоставляющий мощное открытое окружение для инженерных (технических) и научных расчётов.
CeCILL (от «CEA CNRS INRIA Logiciel Libre») — это лицензия на свободное программное обеспечение, адаптированная к интернациональному законодательству и законодательству Франции, подобная GNU General Public License и сохраняющая совместимость с ним.
Система m линейных алгебраических уравнений с n неизвестными (или, линейная система, также употребляется аббревиатура СЛА́У) в линейной алгебре — это система уравнений вида

Обыкнове́нные дифференциа́льные уравне́ния (ОДУ) — это дифференциальное уравнение вида

где — неизвестная функция (возможно, вектор-функция, тогда, как правило, тоже вектор-функция со значениями в пространстве той же размерности; в этом случае говорят о системе дифференциальных уравнений), зависящая от независимой переменной, штрих означает дифференцирование по. Число (порядок старшей производной, входящей в данное уравнение) называется порядком дифференциального уравнения (1)

1. Система компьютерной математики SciLab: История разработки. Возможности и ключевые особенности. Достоинства и недостатки.

История
С 1994 года распространяется вместе с исходным кодом через Интернет. В 2003 году для поддержки Scilab был создан консорциум Scilab Consortium. Сейчас в него входят 25 участников, в том числе Mandriva, INRIA и ENPC (Франция).

Scilab содержит сотни математических функций, и есть возможность добавления новых, написанных на различных языках (C, C++, Fortran и т. д.). Также имеются разнообразные структуры данных (списки, полиномы, рациональные функции, линейные системы), интерпретатор и язык высокого уровня.
Scilab был спроектирован как открытая система, и пользователи могут добавлять в него свои типы данных и операции путём перегрузки.
В системе доступно множество инструментов:
2D и 3D графики, анимация
Линейная алгебра, разреженные матрицы (sparse matrices)
Полиномиальные и рациональные функции
Интерполяция, аппроксимация
Симуляция: решение ОДУ и ДУ
Scicos: гибрид системы моделирования динамических систем и симуляции
Дифференциальные и не дифференциальные оптимизации
Обработка сигналов
Параллельная работа
Статистика
Работа с компьютерной алгеброй
Интерфейс к Fortran, Tcl/Tk, C, C++, Java, LabVIEW
Scilab имеет схожий с MATLAB язык программирования. В состав пакета входит утилита, позволяющая конвертировать документы Matlab в Scilab.
Scilab позволяет работать с элементарными и большим числом специальных функций (Бесселя, Неймана, интегральные функции), имеет мощные средства работы с матрицами, полиномами (в том числе и символьно), производить численные вычисления (например, численное интегрирование) и решение задач линейной алгебры, оптимизации и симуляции, мощные статистические функции, а также средство для построения и работы с графиками.
Для численных расчётов используются библиотеки Lapack, LINPACK, ODEPACK, Atlas и другие.
В состав пакета также входит Scicos — инструмент для редактирования блочных диаграмм и симуляции (аналог simulink в пакете MATLAB). Имеется возможность совместной работы Scilab с программой LabVIEW.

Ключевые особенности
Отличия от некоторых коммерческих программ:
Бесплатность.
Свободность (с версии 5.0).
Маленький размер — дистрибутив 4 версии занимал менее 20 МБ против более чем двухгигабайтного пакета MATLAB. Инсталлятор 5 версии (5.4.0) увеличился в объёме до 108 МБ.
Возможность запуска в консоли без использования графического интерфейса, в том числе в версии под Windows (в UNIX и Windows версиях MatLab-а эта возможность присутствует тоже). Это позволяет производить автоматизированные вычисления, есть пакетный режим.

Достоинства и недостатки
Поиски достоинства и недостатков на просторе интернета ни чего не дали. Так что расскажу о том, что я заметил сам.
Если говорить о достоинствах, тут самым основным для меня является бесплатность данного пакета, по сравнению с той же Mathematic, когда для выполнения лабораторных работы приходилось искать серийник или crack. Далее — это кросплатформеность, т.к. я больше предпочитаю использовать Gentoo, чем Windows. В принципе большенство достоинств описано в пункте Ключевые особенности.
Из недостатков я вижу только 2: Это нет такой визуализации программирования как в Mathematic, а так же система использует прежде всего численные подходы, для вычисления, что может сказаться на точности.

2. Основы работы в SciLab. Пользовательские и системные переменные. Математические выражения. Коментарии

При написании данной статьи я использую версию scilab-5.3.3 под windows.
SciLab 5.3.3
Операционные системы: Windows, Linux, MacOS
Пользовательские и системные переменные
Год/Дата Выпуска: 2011
Версия: 5.3.3
Разработчик: Free Open Source Software for Numerical Computation
Сайт разработчика: www.scilab.org
Разрядность: 32bit+64bit
Лицензия: CeCILL
По характеристикам данных не нашёл, но могу сказать, что данный пакет грузится без глюков на моём ноутбуке
Процессор: Intel Celeron Dual-Core T3300 2.0 ГГц
Оперативная память 2Гб DDR3
Видеоадаптер Intel GMA4500M

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

Есть 2 варианта работы: 1 — это работа в том же командном окне, 2 — открыть SciNotes(что-то вроде блокнота с подсветкой) где можно написать код, который позднее запустить, результат выполнения появится в командном окне.

Для примера я рассмотрю вывод Hellow world.

Командная строка. Используем функцию вывода на дисплей disp()

При работе в SciNotes вы получите что то похожее

для выполнения кода, надо или нажать на стрелочку в право 🙂 (как во многих средах разработки)
или Выполнение->… без отображение команд
в принципе можно использовать и другие методы выполнения, и не использовать вывод на экран
Результатом выполнения будет:

т.к. можно сказать что интерфейс изучен, далее я буду приводить просто код и результат выполнения

SciLab чувствителен к реестру, т.е. А и а — разные переменные переменные.

a=1,A=3
//Каждая операция начинается с новой строки или через запятую
//коментарии можно оставлять после двух символоф слэш
b=3
c=a+b
disp©

Основные операции:
+ сложение
— вычитание
* умножение
/ деление справа, т.е. x/y = xy^(-1)
\ деление слева, т.е. x\y = x^(-1)y
^ возведение в степень, т.е. x^y
** возведение в степень (эквивалентно ^)
’ эрмитово сопряжение (комплексное сопряжение и транспонирование)

3. Основы работы в SciLab. Функции и их типы. Способы объявления пользовательских функций. Использование функций. Файлы-сценарии и их применение для хранения функций.

Элементарные математические функции.

acos acosd acosh acoshm acosm acot acotd acoth
acsc acscd acsch asec asecd asech asin asind
asinh asinhm asinm atan atand atanh atanhm atanm
cos cosd cosh coshm cosm cotd cotg coth
cothm csc cscd csch sec secd sech sin
sinc sind sinh sinhm sinm tan tand tanh
tanhm tanm
exp expm log log10 log1p log2 logm max
maxi min mini modulo pmodulo sign signm sqrt
sqrtm

y = int8(x) 8-битовое число со знаком [-2^7; (2^7)-1] = [-128; 127]
y = uint8(x) 8-битовое число без знака [0; (2^8)-1] = [0; 255]
y = int16(x) 16-битовое число со знаком [-2^15; (2^15)-1] = [-32768; 32767]
y = uint16(x) 16-битовое число без знака [0; (2^16)-1] = [0; 65535]
y = int32(x) 32-битовое число со знаком [-2^31; (2^31)-1] = [-2147483648; 2147483647]
y = uint32(x) 32-битовое число без знака [0; (2^32)-1] = [0; 4294967295]

iconvert преобразование к целочисленному представлению
inttype определение типа целого числа

простейший способ вызова пользовательской функции:

пример пользовательской фунции:

function y = myfunction ( x )
y = 2 * x
endfunction

Сохраняем её. Далее приведён пример вызова данной функции

4. Определение одномерный и многомерных массивов. Основные действия над массивами.

Пример, как задаётся одномерный массив:

квадратные скобки ”[” и ”]” обозначают начало и конец перечисления
элементов матрицы,
запятой ”,” отделяются элементы матрицы, находящиеся в одной строке,
точка с запятой ”;” разделяет строки матрицы.

size определить размер матрицы
matrix изменить размер матрицы
resize_matrix создать новую матрицу заданного размера и скопировать
в нее элементы из исходной матрицы

Операции над матрицами:
Обращение к элементам матрицы
i = 1; 2, а j = 3; 4
для этого возьмём уже готовую матрицу

—>A = testmatrix (» hilb «, 5)
A =
25. — 300. 1050. — 1400. 630.
— 300. 4800. — 18900. 26880. — 12600.
1050. — 18900. 79380. — 117600. 56700.
— 1400. 26880. — 117600. 179200. — 88200.
630. — 12600. 56700. — 88200. 44100.
—>A(1: 2, 3: 4)
ans =
1050. — 1400.
— 18900. 26880

A матрица целиком
A(. ) матрица целиком
A(i:j,k) элементы матрицы в k-ом столбце с i-ой по j-ую строку
A(i,j:k) элементы матрицы в i-ой строке с j-ого по k-ый столбец
A(i,:) i-ая строка матрицы
A(:,j) j-ый столбец матрицы

Генерация единичной матрицы

Операции над матрицами
+ сложение .+ поэлементное сложение
— вычитание .- поэлементное вычитание
* умножение .* поэлементное умножение
/ деление справа ./ поэлементное деление справа
\ деление слева .\ поэлементное деление слева
^ или * возведение в степень :^ поэлементное возведение в степень
’ эрмитово сопряжение (комплексное сопряжение и транспонирование)
.’ транспонирование без сопряжения

пример умножения числа на еденичную матрицу 2 на 2

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

функции работы с матрицами
chol разложение Холесского
companion сопровождающая матрица
cond число обусловленности
det определитель матрицы
inv обратная матрица
linsolve решение систем линейных уравнений
lsq метод наименьших квадратов
lu LU-разложение с выбором опорного элемента
qr QR-разложение
rcond обратное число обусловленности
spec собственные значения и векторы
svd разложение по сингулярным числам матрицы
testmatrix генерация специальных матриц (Гильберта, Франка и др.)
trace след матрицы

6. Определение одномерный и многомерных массивов. Решение СЛАУ. Символьные массивы и операции над ними

Текст файла–сценария с решением задачи по формулам Крамера

—>A=[2 1 -5 1;1 -3 0 -6;0 2 -1 2;1 4 -7 6];//Матрица коэффициентов
—>b=[8;9;-5;0]; //Вектор свободных коэффициентов
—>//Первая вспомогательная матрица
—>A1=A;A1(:,1)=b;
—>//Вторая вспомогательная матрица
—>A2=A;A2(:,2)=b;
—>//Третья вспомогательная матрица
—>A3=A;A3(:,3)=b;
—>//Четвертая вспомогательная матрица
—>A4=A;A4(:,4)=b;
—>//Главный определитель отличен от нуля
—>D=det(A);
—>//Определители вспомогательных матриц
—>d(1)=det(A1);
—>d(2)=det(A2);
—>d(3)=det(A3);
—>d(4)=det(A4);
—>//Вектор неизвестных
—>x=d/D
x =
3.
— 4.
— 1.
1.
—>//Проверка
—>P=A*x-b
P =
0.
0.
— 8.882D-16
2.665D-15

—>A=[2 -1 1;3 2 -5;1 3 -2];
—>b=[0;1;4];
—>//Приведение расширенной матрицы к треугольному виду
—>C=rref([A b]);
—>//Выделение последнего столбца из матрицы,
—>//x — решение системы
—>x=C(1:3,4:4)
x =
0.4642857
1.6785714
0.75
—>A*x //Проверка
ans =
— 5.551D-16
1.
4.

7. Численное интегрирование. Подходы к интегрированию. Интегрирование функций заданных пользователем

Численное интегрирование (историческое название: (численная) квадратура) — вычисление значения определённого интеграла (как правило, приближённое). Под численным интегрированием понимают набор численных методов отыскания значения определённого интеграла.
Интегрирование по методу трапеций
проинтегрируем функцию, корень из 2*x-1 на отрезке от 1 до 10 с разбиением в 1 шаг

—>x=1:10
x =
1. 2. 3. 4. 5. 6. 7. 8. 9. 10.
—>y=sqrt(2*x-1)
y =
column 1 to 6
1. 1.7320508 2.236068 2.6457513 3. 3.3166248
column 7 to 10
3.6055513 3.8729833 4.1231056 4.3588989
—>inttrap(x,y)
ans =
27.211585

Квадратурные формулы Ньютона Котеса

8.Численное дифференцирование. Подходы к дифференцированию.

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

9. Решение ОДУ средствами SciLab. Функции применяемые для решение ОДУ. Решение краевых задач.

Существует 4 способа для решения ОДУ:
1. С помощью команды ode, которая является солвером для решения обыкновенного
дифференциального уравнения.
2. С помощью команды odedc, которая вычисляет решение смешанной дискретно-
непрерывной системы.
3. Команда dassl, которая дает решение неявно выраженного дифференциального уравнения.
4. С помощью команды impl, которая дает решение неявно выраженного линейного
дифференциального уравнения.

—>y0=1;
—>t0=1;
—>t=1:0.01:1.5;
—>deff(«[ydot]=f(t,y)»,«ydot=y^(1/3)*t»)
—>y=ode(y0,t0,t,f);
—>y_exact=((t^2+2)/3)^(1.5);// это функция точного решения для сравнения
—>my_er=y-y_exact;
—>plot(t,y-y_exact) // это график ошибки вычисления от аргумента t

результатом является такой график

за одно можно увидеть Графическое окно. Построение графиков будет подробно рассмотрено далее.

10. Построение двухмерных графиков в системе SciLab. Основные функции и типы графиков.

Функция plot
рассмотрим пример:

Функция plot2d
Рассмотрим функцию опять на примере:

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

Функция polarplot
Служит для построения графика в полярных координатах

получается ромашка

параметры похожи как и в случае с plot2d

11. Построение трёхмерных графиков в системе SciLab. Основные функции и типы графиков.

Существует 4 способа построение графика:
Способ 1.
С помощью команды plot3d. Команда создает 3D график по точкам, заданным матрицами
x, y и z.
Способ 2.
С помощью команды plot3d1. Команда создает 3D график по точкам, заданным
матрицами x, y и z с помощью уровней цвета. Вещь в общем избыточная: величина
координаты z дополнительно еще и покрашена, в зависимости от принимаемого значения
z.
Способ 3.
С помощью команды fplot3d. Это аналог команды fplot3d, но изображаемая поверхность
задана с помощью внешней функции.
Способ 4.
С помощью команды fplot3d1. Это аналог команды plot3d1, но изображаемая поверхность
задана с помощью внешней функции.
Синтаксис этих команд смотри с помощью help.

забыл указать, что в графическом окне есть возможность экспорта данных, т.е. сохранить картинку

Результат такой же как в примере 1

результат такой же, как в примере 2

12. Задача полиномов в SciLab. Символьные операции с полиномами. Решение алгебраический уравнений. Сравнение функций fsolve и roots.

Рассмотрим на примере решения уравнения 2x^4-8x^3+8x^2-1=0

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

Надеюсь данная статья послужит толчком, для дальнейшего изучения SciLab или решения своих задач/лабараторных

Вычисление интегралов, решение уравнений и систем

Краткие теоретические сведения

В Scilab вычисление определенного интеграла методом трапеций реализовано функцией

где x –вектор значений аргумента подынтегральной функции на отрезке интегрирования, y –вектор значений, полученных при вычислении подынтегральной функции для элементов вектора x.

Например, для вычисления нужно выполнить следующий набор команд:

Для вычисления определенного интеграла с использованием алгоритма квадратурных формул предназначена функция

Integrate(fun, x, a, b, ,er1 ,er2),

где fun– подынтегральная функция в символьном виде, x – переменная интегрирования в символьном виде, a, b –пределы интегрирования, er1,er2 –абсолютная и относительная погрешности интегрирования (необязательные параметры).

Например, для вычисления нужно задать:

Универсальная команда интегрирования:

[integral,err]=intg(a, b, name ,er1,er2),

где a, b –пределы интегрирования, name – имя подынтегральной функции (может быть задана с помощью внешней функции), er1,er2 –абсолютная и относительная погрешности интегрирования (необязательные параметры). Функция intgвозвращает значение интеграла (integral) и погрешность вычислений (err).

Внешнюю функцию можно задать командой

deff(‘переменная=имя функции(параметр)’, ‘символьное представление функции’) Например, deff(‘y=F(x)’, ‘y=2*x./(sin(x)+1.5)’).или

function переменная = имя функции(параметр-аргумент функции)

операторы, вычисляющие значение функции

endfunctionНапример, function y=f(t) y=t^2/sqrt(3+sin(t)) endfunctionилиfunction y=f(t),y=t^2/sqrt(3+sin(t)),endfunction

Ниже приводится пример вычисления интеграла .

Для решения нелинейных уравнений в Scilab используется функция

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

Например, для решения уравнения для начального приближения нужно выполнить следующие команды:

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

где а – вектор коэффициентов перед неизвестными полинома размерностью n+1(n – порядок полинома).

Результатом работы этой функции будет вектор корней полинома размерностью n.

Пример решения полиномиального уравнения приведен ниже.

Для уравнения два корня – комплексные.

Функция rootsможет также принимать в качестве параметра полином, созданный функциейpolyи представляющий собой левую часть уравнения :

где a –вектор коэффициентов полинома записанных в обратном порядке, x — символьная переменная, f1– символьная переменная, принимающая значения ‘c’ или ‘r’ (roots или coeff).

Например, чтобы создать полином , нужно использовать команду

Фрагмент программы p=poly([3 -2 1 0 4],’x’,’c’)disp(p)

Результат

Для решения уравнения можно выполнить следующие команды:

Для решения систем линейных уравнений в Scilab есть следующие способы:

— применение операции левого матричного деления;

— использование обратной матрицы.

Если задана система линейных алгебраических уравнений вида:

где А – матрица коэффициентов перед неизвестными системы, В – вектор свободных членов, то решение системы может быть найдено в виде:

То же самое решение может быть получено с помощью обратной матрицы, например:

Например, решить систему уравнений

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

где x0 –вектор начальных приближений для неизвестных, f –функция, определяющая систему

Например, решение системы можно выполнить следующим образом:

Задание 1. Вычисление определенного интеграла

Постановка задачи. Вычислить числовое значение интеграла от этой функции в заданных пределах интегрирования методом трапеций, методом квадратурных формул и с помощью функции intg.

Шаг 1. Создадим вектор X, значения которого будут изменяться от 2,1 до 4,3 с шагом 0.01.

Шаг 2. Создадим вектор Y, каждое значение которого вычисляется по формуле .

Шаг 3. Применим команду inttrap(X, Y).

Шаг 4. Используем функцию integrate,задав подынтегральную функцию в символьном виде.

Шаг 5. Определим внешнюю функцию с помощью команды deffили конструкции function

Шаг 6. Выведем результаты, используя команду disp.

Индивидуальные задания приведены в таблице 3.4.

Задание 2. Поиск корней уравнения, графическая интерпретация

Постановка задачи. Найти корень уравнения для заданного начального приближения. Выполнить графическую интерпретацию результата.

Шаг 1. Определим внешнюю функцию с помощью команды deffили конструкции function

Шаг 2. Найдем корень уравнения с помощью функции fsolve,подставив в качестве первого параметра заданное начальное приближение.

Шаг 3. Выведем результат, используя команду disp.

Шаг 4. Выполним графическую интерпретацию результата. Для этого зададим аргумент функции из левой части уравнения таким образом, чтобы найденный корень попадал в диапазон между первым и последним элементом вектора. Построим график функции из левой части уравнения с помощью plot.Построим также линию y=0 и отметим точку с абсциссой, равной корню, и ординатой, равной значению функции для корня.

Программа	Результат выполнения
deff(‘y=F(x)’, …’y=sin(x)-cos(x.^2)-sin(2*x)’)root=fsolve(1,F)disp(root)x=0.5:0.01:2.5plot(x,F(x),’-b’,root,F(root),’xr’,x,0,’-k’)	1.1695683

Индивидуальные задания приведены в таблице 3.5

Вариант	Уравнение	Начальноеприближение
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.

Задание 3. Поиск корней полиномиального уравнения, графическая интерпретация

Постановка задачи. Найти все корни полиномиального уравнения. Выполнить графическую интерпретацию для одного из найденных действительных корней.

Шаг 1. Создадим вектор коэффициентов полинома в левой части уравнения (или полином с помощью poly)

Шаг 2. Найдем корни уравнения с помощью функции roots.

Шаг 3. Выведем результат, используя команду disp.

Шаг 4. Выполним графическую интерпретацию результата. Для этого зададим аргумент функции из левой части уравнения таким образом, чтобы выбранный действительный корень попадал в диапазон между первым и последним элементом вектора. Построим график функции из левой части уравнения с помощью plot.Построим также линию y=0 и отметим точку с абсциссой, равной корню, и ординатой, равной значению функции для корня.

Программа	Результат выполнения
v=[2 0 4 -6 -3]R=roots(v)disp(R)root=R(3)x=0.5:0.01:2y=2*x.^4+4*x.^2-6*x-3F_root=2*root^4+4*root^2-6*root-3plot(x,y,’-b’,root,F_root,’xr’,x,0,’-k’)	— 0.4129576 + 1.7282075i — 0.4129576 — 1.7282075i 1.2164706 — 0.3905555

Индивидуальные задания приведены в таблице 3.6.

Вариант	Уравнение
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.

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

Постановка задачи. Решить систему линейных уравнений.

Шаг 1. Создадим матрицу коэффициентов при неизвестных

Шаг 2. Создадим вектор свободных членов.

Шаг 3. Умножим матрицу, обратную к матрице коэффициентов, на вектор свободных членов (или применим операцию левого матричного деления).

Шаг 4. Выведем результат, используя команду disp.

Индивидуальные задания приведены в таблице 3.7.

№	Система уравнений	№	Система уравнений
1.		2.
3.		4.
5.		6.
7.		8.
9.		10.
11.		12.
13.		14.
15.		16.

Статьи к прочтению:

Метод Крамера за 3 минуты. Решение системы линейных уравнений — bezbotvy

Похожие статьи:

Приближенное вычисление интегралов Приближённое вычисление определённого интеграла основано на геометрическом смысле интеграла и сводится к приближённому…

Задание: графически и численно решить систему нелинейных алгебраических уравнений, на примере поиска точек пересечения двух функций. Исходные данные:…

Решение дифференциальных уравнений в Scilab

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

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

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

Дифференциальные уравнения n-го порядка

Дифференциальным уравнением n-го порядка называется соотношение вида

Решением дифференциального уравнения является функция \( x(t) \), которая обращает это уравнение в тождество. Дифференциальные уравнения имеют бесконечное множество решений.

Чтобы однозначно определить решение некоторой задачи, необходимо задавать начальные условия. Количество таких условий должно совпадать с порядком уравнения. А решение, найденное для таких условий называется решением задачи Коши.

6.2 Пример решения оду в Scilab

Для решения обыкновенных дифференциальных уравнений (ОДУ) в Scilab используется
функция y = ode([type], y0, t0, t, func).
Разберём, что обозначает каждый из параметров у этой функции.

type — необязательный строковый параметр, с помощью которого можно выбрать метод решения ОДУ. Обычно этот параметр опускается.
t0 — скаляр начальный момент отрезка интегрирования. Обычно \( t0 = 0 \).
y0 — начальные условия. Отметим, что \(y0 \) это вектор, размерность которого совпадает с порядком ОДУ. Так, для ОДУ 2-го порядка необходимо задать значения в начальный момент времени для функции и её производной, т.е. использовать запись \(y0 = [0.1, 0.3]\).
t — вектор, задающий узлы сетки, в которых ищем решение. Как правило, вектор t задается следующим образом t=t0:d:tmax, где \(t0\) — начальный момент отрезка интегрирования, \(d\) — шаг дискретизации, \(tmax\) — конечный момент отрезка интегрирования.
func — пользовательская функция, определяющая правую часть уравнения.
y — вектор решений.

Рассмотрим использование функции ode() на примере решения следующей задачи.

Найти угловую скорость \( \omega(t) \) твердого тела вокруг неподвижной оси, если заданы начальная угловая скорость тела \( \omega_0=1 \) и угловое ускорение \( \varepsilon(t) = 2+0.5sin(t) \). Найти значение \( \omega(t) \) в момент времени \( t_1=1.3 \). Построить график функций.

Рисунок 9. Сравнение графиков угловой скорости точки, по функции, найденной аналитическим и численным интегрированием. Красным кружком обозначено значение скорости точки в момент времени t=1.3c.

Угловая скорость точки может быть найдена из дифференциального уравнения

В нашем случае по условию задачи указаны следующие параметры для функции ode():
\(t0=0 \) начальный момент времени,
\(y0 = 1\) начальное условие одно, т.к. порядок уравнения равен 1, это начальная угловая скорость тела \( \omega_0=1 \)
t=t0:h:tmax — вектор, задающий узлы сетки.

Результат работы программы представлен на рис.9, исходный код на листинге 15.

Отметим, что в данной программе фигурирует пользовательская функция aomega(t), представляющая собой угловую скорость, найденную аналитически с указанным начальным условием.

Функция \( f = 2t — 0,5cos(t) + 1,5 \) единственное решение дифференциального уравнения, удовлетворяющее заданным начальным условиям, то есть, функция \(f\) решение задачи Коши.

6.3 Cистемы дифференциальных уравнений

Для решения систем ОДУ в Scilab используется та же функция y = ode([type], y0, t0, t, func). Однако важным требованием является запись исследуемой системы в нормальной форме Коши.

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

\begin x_1′ = f_1(t,x_1. x_n) \\ x_2′ = f_2(t,x_1. x_n) \\ . \\ x_n’ = f_n(t,x_1. x_n) \end

Решением системы ОДУ является вектор x(t), который обращает это уравнение в тождество. Размерность вектора равна количеству уравнений в системе.

Стоит отметить, что дифференциальное уравнение n-ой степени может быть представлено в виде системы из n-уравнений первой степени, что позволяет решать задачу Коши для полученной системы ОДУ.

6.4 Примеры поиска решения систем ОДУ в Scilab

Решение системы ОДУ

Рассмотрим решение задачи Коши для системы уравнений

\begin x’ = cos(xy) \\ y’ = sin(x+ty) \end

на интервале \( [0; 10]\) и с начальными условиями \(x(0)=0, y(0)=0 \).

Для поиска решения данной задачи, нам необходимо привести исходную систему к нормальному виду Коши. Для этого введём новые переменные \( (z1, z2) \) и сделаем необходимые переобозначения в исходной системе:

Рассмотрим решение задачи Коши для системы уравнений

\begin z_1 = x \\ z_2 = y \end \begin z_1′ = cos(z_1 z_2) \\ z_2′ = sin(z_1 + t z_2) \end

Код программы, реализующей поиск решения системы дифференциальных уравнений представлен на листинге 16. Обратите внимание, как происходит обращение к компонентам вектора решения \( z(t) \).

Для наглядной реализации сформировано графическое решение на рис. 10.

Рисунок 10. Графическое решение задачи Коши с помощью функции ode().

Решение системы ОДУ в матричной форме

Рассмотрим решение задачи Коши для системы уравнений, заданных в матричном виде

на интервале \( [0; 10] \) и начальными условиями \( y_1(0)=1, y_2(0)=0 \).

Код программы, реализующей поиск решения системы дифференциальных уравнений представлен на листинге 17. Обратите внимание, как происходит обращение к компонентам вектора решения \(y(t)\).

Для наглядной реализации сформировано графическое решение на рис. 11.

Рисунок 11. Графическое решение задачи Коши с помощью функции ode().

Решение ДУ 2-го порядка путём сведения к системе уравнений

Продолжим знакомиться с возможностями функции ode() на примере решения задачи механики на второй закон Ньютона.

Груз находится на пружине жёсткости \( c=12 \) H/м, масса груза \( m=68.7 \) кг. Определить закон движения груза, если на него действует сила \( F=100.5sin(2t) + <2\pi \over 3>\) H.

По 2-му закону Ньютона, движение груза описывается с помощью дифференциального уравнения 2-й степени, которое имеет вид:

и начальными условиями \( x(0)= 0, ;\ \dot(0)= 0 \).

Как известно, ОДУ второго порядка сводится к системе в нормальной форме Коши, состоящей из двух уравнений первой степени. Введём новые переменные \( (z_1, z_2) \) и сделаем необходимые переобозначения в исходной системе:

\begin z_1 = x \\ z_2 = \dot \end \begin \dot_1 = z_2 \\ \dot_2 = <1 \over m>(-cz_1 + F(t)) \end

Код программы, реализующей поиск решения системы дифференциальных уравнений представлен на листинге 18. Обратите внимание, как происходит обращение к компонентам вектора решения \( y(t) \).

Для наглядной реализации сформировано графическое решение на рис. 12. Компонента \( z_1 \) представляет собой координату груза, а компонента \( z_2 \) — скорость груза.

Рисунок 12. Графики движения груза на пружине для компонент \( z_1 и z_2 \) соответственно.