Решение уравнений в SciLab
СОДЕРЖАНИЕ
ВВЕДЕНИЕ
Создание математических моделей различных физических процессов и построение методов решения физических задач является предметом математической физики. Постановка задач математической физики заключается в построении математических моделей, описывающих основные закономерности изучаемых физических явлений. Это моделирование состоит в выводе уравнений, которым удовлетворяют величины, характеризующие данный процесс. При этом исходят из основных физических законов, которые учитывают наиболее существенные черты явления, отвлекаясь от второстепенных характеристик.
Многие задачи механики и
физики приводят к исследованию дифференциальных
уравнений с частными производными.
В ряде случаев приходится иметь
дело с интегральными и интегро-
Автоматизация научных исследований позволяет получать более точные и полные модели исследуемых объектов и явлений, ускорять ход научных исследований и снижать их трудоемкость, изучать сложные объекты и процессы, исследование которых традиционными методами затруднительно, позволяет.
Автоматизированные системы научных исследований (АСНИ) отличаются от других типов автоматизированных систем (АСУ, АСУТП, САПР и т.д.) характером информации, получаемой на выходе системы. Прежде всего, это обработанные или обобщенные экспериментальные данные, но главное - полученные на основе этих данных математические модели исследуемых объектов, явлений или процессов. Адекватность и точность таких моделей обеспечивается всем комплексом методических, программных и других средств системы. В АСНИ могут использоваться также и готовые математические модели для изучения поведения тех или иных объектов и процессов, а также для уточнения самих этих моделей. АСНИ поэтому являются системами для получения, корректировки или исследования моделей, используемых затем в других типах автоматизированных систем для управления, прогнозирования или проектирования.
Существует немало АСНИ, которые используются в различных отраслях научной деятельности. Среди них можно выделить пакет прикладных программ SCILAB, который распространяется свободно.
1. ХАРАКТЕРИСТИКА СИСТЕМЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ SCILAB
Scilab — пакет прикладных математических программ, предоставляющий мощное открытое окружение для инженерных (технических) и научных расчётов.
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.
Программа доступна для различных операционных систем, включая Linux и Microsoft Windows. Возможности Scilab могут быть расширены внешними программами и модулями, написанными на разных языках программирования. Программа имеет открытый исходный код, что позволяет как свободное коммерческое использование и распространение неизменённых версий, так и некоммерческое распространение измененных версий, которые должны включать в себя исходный код. Для коммерческого распространения измёненных версий необходимо согласование с INRIA.
Начиная с версии 5.0 программа распространяется под совместимой с GNU GPL 2 лицензией CeCILL.
Отличия от некоторых коммерческих программ:
- Бесплатность.
- Свободность (с версии 5.0).
- Маленький размер — дистрибутив 4 версии занимал менее 20 МБ против более чем двухгигабайтного пакета MATLAB. Инсталлятор 5 версии (5.3.2) увеличился в объёме до 121 МБ.
- Возможность запуска в консоли без использования графического интерфейса, в том числе в версии под Windows (в UNIX и Windows версиях MatLab-а эта возможность присутствует тоже). Это позволяет производить автоматизированные вычисления, есть пакетный режим.
2. РЕШЕНИЕ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ
2.1.УРАВНЕНИЕ ТЕПЛОПРОВОДНОСТИ
Пусть физическое тело, температура которого в точке с координатами (x, y, z) в момент времени t определяется функцией u(x, y, z, t). Если температура различных частей тела различна, то в теле будет происходить перенос тепла от более нагретых участков тела к менее нагретым. Рассмотрим какую-нибудь поверхность S внутри тела и обозначим ее малый элемент ∆s. Тогда количество тепла, проходящего через элемент поверхности ∆s за время ∆t, будет определяться выражением
(2.1)
где k > 0 – коэффициент теплопроводности, а – нормаль к элементу поверхности ∆s в направлении распространения тепла. Здесь и далее , a ▽u означает градиент функции u. Будем считать, что тело изотропно в отношении теплопроводности, т. е. коэффициент теплопроводности не зависит от направления нормали, а зависит только от координат точки. Кроме того, обозначим поток тепла через единичную площадку через q т. е.
(2.2)
Тогда для произвольного объема V , ограниченного замкнутой поверхностью S, изменение количества тепла за промежуток времени (t1, t2) будет определяться интегралом
(2.3)
где – внутренняя нормаль к поверхности S.
С другой стороны, количество тепла, необходимое для изменения температуры выделенного объема тела на ∆u = u(x, y, z, t2)−u(x, y, z, t1),
равно
(2.4)
где ρ и γ – плотность и теплоемкость данного физического тела.
Если внутри рассматриваемого тела имеются источники тепла плотности F(x, y, z, t), то количество тепла, выделяемого (или поглощаемого) ими в объеме V за промежуток времени (t1, t2), будет равно
(2.5)
Очевидно, что баланс тепла для объема V определяется соотношением Q2 = Q1 + Q3, т. е.
Применяя к поверхностному
интегралу формулу Гаусса-
Получим
(2.7)
где
Так как подынтегральная функция непрерывна, а объем V и промежуток времени (t1, t2) произвольны, то для любой точки (x, y, z) и любого момента времени справедливо соотношение
(2.8)
Это уравнение называется уравнением теплопроводности неоднородного изотропного тела. Если тело однородно, то коэффициенты γ, ρ, k являются постоянными и уравнение теплопроводности принимает вид
(2.9)
где
Если кроме того отсутствуют внутренние источники тепла, то получим однородное уравнение теплопроводности в виде
(2.10)
Для ряда частных задач уравнение теплопроводности может быть сведено к двумерному или одномерному случаю. Так, для очень тонкой однородной пластины уравнение (1.10) принимает вид
(2.11)
В случае тонкого однородного стержня уравнение (2.10) сводится к одномерному уравнению
(2.12)
Следует заметить, что во всех описанных выше случаях тепловой обмен с окружающей средой не учитывается.
Для корректной постановки задачи о нахождении температуры в любой точке тела в каждый момент времени следует к уравнению теплопроводности, как и в случае струны, добавить начальные и граничные условия, которые будут определять начальное распределение температуры внутри тела и тепловой режим на границе. Начальные условия будут иметь вид
(2.13)
Граничные условия могут задаваться несколькими различными способами, например:
1. В каждой точке поверхности M(x, y, z, t) ∈ S задается температура
(2.14)
2. На поверхности тела задается тепловой поток , который определяет граничные условия вида
(2.15)
Таким образом, в общем случае задача о распространении тепла в изотропном теле корректно ставится следующим образом: найти решение уравнения (2.8) , удовлетворяющее начальному условию (2.13) и одному из граничных условий (2.14) или (2.15).
2.2.ИСПОЛЬЗОВАНИЕ МЕТОДА СЕТОК ДЛЯ РЕШЕНИЯ ПАРАБОЛИЧЕСКИХ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ
Пусть дано дифференциальное уравнение (2.12).
Требуется найти функцию в области с границей при заданных краевых условиях. Согласно методу сеток в плоской области строится сеточная область , состоящая из одинаковых ячеек. При этом область должна как можно лучше приближать область . Сеточная область (то есть сетка) состоит из изолированных точек, которые называются узлами сетки. Число узлов будет характеризоваться основными размерами сетки : чем меньше , тем больше узлов содержит сетка. Узел сетки называется внутренним, если он принадлежит области , а все соседние узлы принадлежат сетке . В противном случае он называется граничным. Совокупность граничных узлов образует границу сеточной области .
Сетка может состоять из
клеток разной конфигурации: квадратных,
прямоугольных, треугольных и других.
После построения сетки исходное
дифференциальное уравнение заменяется
разностным уравнением во всех внутренних
узлах сетки. Затем на основании
граничных условий
Замена дифференциального уравнения разностным может быть осуществлена разными способами.
Один из способов аппроксимации состоит в том, что производные, входящие в дифференциальное уравнение, заменяются линейными комбинациями значений функции в узлах сетки по тем или иным формулам численного дифференцирования.
Различные формулы численного
дифференцирования имеют разную
точность, поэтому от выбора формул
аппроксимации зависит качество
аппроксимации
Рассмотрим неоднородное уравнение теплопроводности
, (2.16)
- известная функция. Будем
Заметим, что эту полуполосу всегда можно привести к полуполосе, когда . Уравнение (2.16) будем решать с начальными условиями:
,
- известная функция, и краевыми условиями:
(2.17)
где - известные функции переменной .
Для решения задачи область покроем сеткой .
Узлы сетки, лежащие на прямых , и будут граничными. Все остальные узлы будут внутренними. Для каждого внутреннего узла дифференциальное уравнения (1.16) заменим разностным. При этом для производной воспользуемся следующей формулой:
.
Для производной запишем следующие формулы:
,
,
.
Можем получить три вида разностных уравнений:
, (2.18)
, (2.19)
, (2.20)
.
Разностные уравнения (2.18) аппроксимируют уравнение (2.16) с погрешностью , уравнение (2.19) - с такой же погрешностью, а уравнение (1.7) уже аппроксимирует уравнение (2.16) с погрешностью .
В разностной схеме (2.17) задействованы 4 узла. Конфигурация схемы (2.17) имеет вид:
В схеме (2.19) также участвуют 4 узла, и эта схема имеет вид:
В схеме (2.20) участвуют 5 узлов, и эта схема имеет вид:
Первая и третья схемы - явные, вторая схема неявная. В случае явных схем значения функции в узле очередного слоя можно найти, зная значения в узлах предыдущих слоев. В случае неявных схем для нахождения значений решения в узлах очередного слоя приходится решать систему уравнений.
Для узлов начального (нулевого) слоя значения решения выписываются с помощью начального условия:
Для граничных узлов, лежащих на прямых и , заменив производные по формулам численного дифференцирования, получаем из граничных условий следующие уравнения:
(2.21)
Уравнения (2.21) аппроксимируют граничные условия с погрешностью , так как используем односторонние формулы численного дифференцирования. Погрешность аппроксимации можно понизить, если использовать более точные односторонние (с тремя узлами) формулы численного дифференцирования.
Присоединяя к системе разностных уравнений, записанных для внутренних узлов, начальные и граничные условия для разностной задачи получим полные разностные схемы трех видов. Для проведения вычислений самой простой схемой оказывается первая: достаточно на основании начального условия найти значения функции в узлах слоя , чтобы в дальнейшем последовательно определять значения решения в узлах слоев и т.д.
Третья схема также весьма проста для проведения вычислений, но при ее использовании необходимо кроме значений решения в узлах слоя найти каким-то образом значения функции и в слое . Далее вычислительный процесс легко организовывается. В случае второй схемы, которая является неявной, обязательно приходится решать систему уравнений для нахождения решения сеточной задачи.
С точки зрения точечной аппроксимации третья схема самая точная.
Введем в рассмотрение параметр . Тогда наши разностные схемы можно переписать, вводя указанный параметр. При этом самый простой их вид будет при .
В любом случае согласно методу сеток будем иметь столько уравнений, сколько имеется неизвестных (значения искомой функции в узлах). Число неизвестных равно числу всех узлов сетки. Решая систему уравнений, получаем решение поставленной задачи.
Разрешимость этой системы для явных схем вопросов не вызывает, так как все действия выполняются в явно определенной последовательности. В случае неявных схем разрешимость системы следует исследовать в каждом конкретном случае. Важным вопросом является вопрос о том, на сколько найденные решения хорошо (адекватно) отражают точные решения, и можно ли неограниченно сгущая сетку (уменьшая шаг по осям) получить приближенные решения, сколь угодно близкие к точным решениям? Это вопрос о сходимости метода сеток.
На практике следует применять сходящиеся разностные схемы, причем только те из них, которые являются устойчивыми, то есть при использовании которых небольшие ошибки в начальных или промежуточных результатах не приводят к большим отклонениям от точного решения. Всегда следует использовать устойчивые разностные схемы, проводя соответствующие исследования на устойчивость.
Первая из построенных выше разностных схем в случае первой краевой задачи будет устойчивой при . Вторая схема устойчива при всех значениях величины . Третья схема неустойчива для любых , что сводит на нет все ее преимущества и делает невозможной к применению на ЭВМ.
Явные схемы просты для
организации вычислительного
Задача и листингн программы
Найти распределение температуры в стержне длиной L, начальная температура
стержня задается произвольной функцией f(x). Температуры концов стержня равны
u0,t=U1=const , uL,t=U2=const. На боковой поверхности стержня
происходит теплообмен по закону Ньютона со средой температура которой равна ,
u0 .
//Начальное условие
function y=fi(x)
y=exp(0.15*x)
endfunction
function [u,x,t,r,k]=neiav(N,K,L,T,a,h,
//Функция решения
// помощью неявной разностной схемы.
//N - количество участков, на которые разбивается интервал по13
х
// (0,L); K - количество участков, на которые разбивается
// интервал по t (0,T);
//a, h - параметры диффренциального уравнения
теплопроводности;
// eps - точность решения СЛАУ) методом Зейделя;
// U1 - температура на левом конце стержня;
// U2 - температура на правом конце стержня;
// Функция neiav возвращает:
// u - матрицу решений в узлах сетки, массив x, массив t,
// r - точность решения системы методом Зейделя,
// к - количество итераций при решении системы методом
Зейделя.
// Вычисляем шаг по х
hx=L/N;
// Вычисляем шаг по t
delta=T/K;
// Формируем массив x и первый столбец матрицы решений U
// из начального условия
for i=1:N+1
x(i)=(i-1)*hx;
u(i,1)=fi(x(i));
end
// Формируем массив t, первую и последнюю строку матрицы
// решений U из граничных условий
for j=1:K+1
t(j)=(j-1)*delta;
u(1,j)=U1;
u(N+1,j)=U2;
end
// Определяем матрицу ошибок R и заполняем ее нулями
R(N+1,K+1)=0;
// Вычисляем коэффициент gamma
gam=a^2*delta/hx^2;
r=1;
k=0;
//Цикл while для организации итерационного процесса при
решении
//системы уравнений методом Зейделя с точностью eps
while r>eps
// Вычисление матрицы ошибок R во внутренних точках
//и пересчет значений u во внутренних точках при решениии СЛАУ
// (12.17) методом Зейделя
for j=2:K+1
for i=2:N
V=gam*(u(i-1,j)+u(i+1,j))/(1+
(1+2*gam+delta*hx)+delta*h*u0/
R(i,j)=abs(V-u(i,j));
u(i,j)=V;
end14
end
// Поиск максимума в матрице ошибок
r=R(1,1);
for i=1:N+1
for j=1:K+1
if R(i,j)>r
r=R(i,j);
end
end
end
// Увеличение количества итераций
k=k+1;
end
endfunction
//Вызов функции решения задачи
[U,X,T,R,K]=neiav(50,200,5,3,
//Построение графика функции
surf(X,T,U');
title('Example 12.2');
xlabel('X');
ylabel('T');
График решения задачи
Заключение
В своей дипломной работе я бы хотел подробно изучить и расмотреть задачу решения уравнения теплорроводимости в программе Scilab.
Многим из тех, кто сталкивается с научными и инженерными расчётами часто приходится оперировать наборами значений, полученных опытным путём или методом случайной выборки. Как правило, на основании этих наборов требуется построить функцию, на которую могли бы с высокой точностью попадать другие получаемые значения.
В ходе проделанной работы была изучена литература, описывающая численные методы, программирование в среде Free Pascal и SciLab.
Также изучена
литература, описывающая физическую
интерпретацию данного
С помощью среды SciLab были проверены вычисления, проводимые в прикладной программе. Эта проверка показала достаточно высокую точность вычислений прикладной программы. Значения, полученные в среде SciLab, отличаются от вычислений прикладной программы примерно на одну тысячную. Это позволяет судить о достаточно высокой точности вычислений прикладной программы и соответственно о точности вычислений интеграла методом трапеций, который использовался в прикладной программе.