Обзор численных методов в математическом моделировании

 

СОДЕРЖАНИЕ

Введение…………………………………………….………….……….............   

1. Обзор численных методов в математическом   моделировании…………………………………………………..............
1. Анализ численных методов в математическом моделировании…………
2. Аппроксимация и интерполяция в MathCAD………………………………
3. Численные методы решения ОДУ и их реализации в MathCAD…………
2. АЛГОРИТМИЧЕСКИЙ АНАЛИЗ   ЗАДАЧИ……………………......…….
1. Полная постановка задачи…………………………………………...……...
2. Описание математической модели…………………………………...…….
3. Анализ исходных и результирующих данных……………………...……..
4. Графическая схема алгоритма и её описание………………………..........
3. РЕАЛИЗАЦИЯ МОДЕЛИ В MathCad………………………………..…….
1. Описание реализация базовой модели……………………………….........
2. Описание исследований……………………………………………..……..
3. Выводы по результату исследований………………………………..……

ЗАКЛЮЧЕНИЕ……………………………………………………………..........

СПИСОК ИСПОЛЬЗОВАННОЙ ЛИТЕРАТУРЫ………………….……….

ПРИЛОЖЕНИЕ А. Базовая модель……………………………………….......

ПРИЛОЖЕНИЕ В. Исследование……………………………………………..

ПРИЛОЖЕНИЕ С. Индивидуальное задание……………………………….

 

 

 

 

 

 

 

 

 

 

 

         

 

 

1. обзор численных методов в математическом моделировании

 

1.1 Анализ численных методов в математическом моделировании

 

Наиболее распространенными из численных методов, применяемых в математическом  моделировании, являются метод Эйлера и метод Рунге-Кутта.

Метод Эйлера

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

Пусть требуется решить задачу Коши: найти решение дифференциального уравнения

у'=(х, у),                                               (1.1.1)  

 удовлетворяющее условию

у=у0    при    х=х0, т. е. у(х0) =у0.                       (1.1.2)  

При численном решении уравнения (1.1.1) задача ставится так: точках х0, x1, х2, ..., хп найти приближения уп для значений точного решения у(хп). Разность x=xn+1—xn = h называется шагом сетки. Во многих случаях принимают величину h постоянной, тогда

             Хn=Х0=nh(n=0,1,2,…)                                 (1.1.3)

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

 

 

Метод Эйлера состоит в непосредственной замене производной разностным отношением по приближенной формуле:

         (1.1.4)

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

 

         (1.1.5)  

            В частности, при x=x1 получаем y1=y0+ht(x0, y0). Равенство (1.1.4) означает, что на отрезке [х0, xo+h] искомую интегральную кривую у=у(х) приближенно заменяют прямолинейным отрезком, выходящим из начальной точки М0(х0, у0) с угловым коэффициентом f(xQ, у0). Аналогично находим приближенное значение y2: y2 = y1+hf(x1,y1).

Для точки xn = xo+nh получаем

                            (1.1.6)

 Вычисление приближений уп искомого решения у(х) по формуле (1.1.4) представляет собой обыкновенный метод Эйлера. Этот метод дает весьма грубое приближение решения задачи Коши. Он обычно используется в случае, когда необходимо получить примерное представление о решении на небольшом промежутке.

Погрешность обыкновенного метода Эйлера оценивается формулой

                     (1.1.7)

Метод Рунге – Кутта

Метод Рунге-Кутта является одним из наиболее употребительных численных методов повышенной точности. Низкая точность метода Эйлера связана в первую очередь с тем, что остаточный член формулы Эйлера велик. Очевидно, что для уменьшения погрешности вычисления необходимо увеличить количество учитываемых членов в формуле Тейлора. Наиболее распространенным является метод Рунге-Кутта 4-го порядка, в котором учтены производные до 4-го порядка включительно. Метод Эйлера можно рассматривать как метод Рунге-Кутта 1-го порядка. Метод Рунге-Кутта требует большого объёма вычислений, однако расчёт оказывается более точным, чем расчёт по методу Эйлера с тем же шагом.

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

Методы Рунге-Кутта обладают следующими свойствами:

1. Эти методы являются одноступенчатыми: чтобы найти уm+1, нужна информация  о предыдущей точке xm , ym.

2. Они согласуются с рядом Тейлора вплоть до членов порядка hp, где степень р    различна для различных методов и называется порядковым номером или    порядком метода.

3. Они не требуют вычисления производных от f (x,y), а требуют вычисления   самой  функции.

Наиболее точный метод решения – метод Рунге-Кутта 4 порядка, один из самых употребляемых методов интегрирования дифференциальных уравнений.

Этот  метод:

- является одноступенчатым  и одношаговым;

 

- требует  информацию  только  об одной  точке;

- имеет небольшую погрешность;

- значение функции  рассчитывается  при каждом шаге.

Формулы, описывающие классический метод Рунге-Кутта четвертого порядка, состоят из следующих пяти соотношений :

ym+1=ym+h/6(R1+2R2+2R3+R4) ,    (1.1.8)

где

R1=f(xm,ym)                                       (1.1.9)

R2=f(xm+h/2,ym+hR1/2)      (1.1.10)

R3=f(xm+h/2,ym+hR2/2)     (1.1.11)

R4=f(xm+h/2,ym+hR3/2)      (1.1.12)

Ошибка ограничения для этого метода равна et=kh5  Значит формулы

 1.1.8 – 1.1.12 описывают метод Рунге-Кутта четвертого порядка. Однако при использовании этого метода функцию необходимо вычислять четыре раза.

 

2. Апроксимация и интерполяция в MathCAD

 

Интерполяция использует значения некоторой функции, заданные в ряде точек, чтобы предсказать значения функции между ними. В MathCAD можно соединять точки данных прямыми линиями (линейная интерполяция) или соединять их отрезками кубического полинома (кубическая сплайн-интерполяция).

Функции интерполяции определяют кривую, точно проходящую через заданные точки. Из-за этого результат очень чувствителен к ошибкам данных. Кроме того, каждый элемент массива, который используется в любой из функций, описанных в этом разделе, содержит определенное значение. Поскольку MathCAD присваивает значение 0 любым элементам, которые явно не определены. 

Для построения интерполяции в MathCAD имеются несколько встроенных функций, позволяющих "соединить" точки выборки данных (xi,yi) кривой разной степени гладкости. По определению, интерполяция означает построение функции D(х), аппроксимирующей зависимость у(х) в промежуточных точках. Поэтому интерполяцию еще по-другому называют аппроксимацией. В точках xi значения интерполяционной функции должны совпадать с исходными данными, т. е. A(xi)=y(xi).

Самый простой вид интерполяции - линейная, которая представляет искомую зависимость А(х) в виде ломаной линии. Интерполирующая функция А(х) состоит из отрезков прямых.

В MathCAD для построения линейной интерполяции служит встроенная функция

linterp(х, у, t),

где

х - вектор действительных данных аргумента;

у - вектор действительных данных значений того же размера;

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

Элементы вектора х должны быть определены в порядке возрастания, т. е. Х1<Х2<Х3<. . . <XN.

Для проведения кубической сплайн-интерполяции MathCAD предлагается три встроенные функции (VX, VY – вектора узловых точек):

-cspline(VX, VY) – возвращает вектор вторых производных (VK) при приближении в опорных точках к кубическому полиному;

-pspline(VX, VY) – возвращает вектор вторых производных (VK) при приближении в опорных точках к параболической кривой;

-lspline(VX, VY) – возвращает вектор вторых производных (VK) при приближении в опорных точках к прямой.

Интерполирующая функция строится с помощью стандартной функции interp, имеющей следующий общий вид:

interp(VK,VX, VY, x),

где

VK – вектор вторых производных  сплайна в опорных точках;

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

MathCAD позволяет проводить линейную регрессию общего вида, в которой аппроксимирующая функция задается линейной комбинацией функций, причем сами функции fi(x) могут быть нелинейными:

                                                                                                                 (1.2.1)

Линейная регрессия общего вида реализуется с помощью функции linfit:

linfit(VX, VY, F),

где

VX, VY – координаты исходных точек;

F  - вектор, содержащий функции fi(x) , записанные в символьном виде.

Функция linfit еще называется функцией аппроксимации по методу наименьших квадратов.

Результатом работы функции linfit является вектор коэффициентов К, при котором среднеквадратичная погрешность приближения  исходных точек с координатами VX, VY, минимальна.

Вектор VX должен быть возрастающим.

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

Для этого используется встроенная функция genfit, имеющая следующий общий вид:

genfit(VX,VY,VS,F),

где

VS - вектор, который задает начальные приближения элементов вектора K, рассчитываемых итерационным способом;

F - вектор, который содержит искомую функцию и ее частные производные по параметрам Ki в аналитическом виде:

                                                                                                          

                                                                                                                  (1.2.2)

 

 

3. Численные методы решения ОДУ и их реализации в MathCAD

 

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

- rkfixed – функция для решения ОДУ и систем ОДУ методом Рунге–Кутта четвертого порядка с постоянным шагом;

- Rkadapt – функция решения ОДУ и систем ОДУ методом Рунге–Кутта с переменным шагом;

- Odesolve – функция, решающая ОДУ блочным методом.

Ниже приведено описание стандартной функции rkfixed с указанием параметров функции:

rkfixed(y, x1, x2, p, D)

где

y – вектор начальных  условий из k элементов (k – количество  уравнений в системе);

x1 и x2 – левая и правая  границы интервала, на котором  ищется решение ОДУ или системы  ОДУ;

p – число точек внутри  интервала (x1, x2), в которых ищется  решение;

D – вектор, состоящий  из k-элементов, который содержит  первую производную искомой функции  или первые производные искомых  функций, если речь идет о решении  системы.

Результатом работы функции  является  матрица из p+1 строк, первый столбец которой содержит точки, в которых получено решение, а остальные столбцы – сами решения.

При решении дифференциального уравнения первого порядка нужно создать вектор начальных условий из одного элемента Y1, который затем используется при формировании вектора-функции правой части дифференциального уравнения. При обращении к функции rkfixed указывается имя вектора Y, границы интервала, на котором ищется решение уравнения, например, (0;5), количество точек, в которых ищется решение – 100, вектор-функция, описывающая правую часть дифференциального уравнения – D. В результате получается матрица, в первом столбце которой содержатся значения аргумента искомой функции, во втором – значения самой результирующей функции. При построении графика функции первый столбец полученной матрицы указывается как аргумент, второй столбец – как функция.

 

При решении системы дифференциальных уравнений нужно создать вектор начальных условий из двух элементов, например, вектор v, который затем используется при формировании вектора-функции правой части дифференциального уравнения. При обращении к функции rkfixed указывается имя вектора v, и границы интервала, на котором ищется решение уравнения, например, (0;5), количество точек, в которых ищется решение – 100, вектор-функция, описывающая правую часть дифференциального уравнения – D. В результате получается матрица, в первом столбце которой содержатся значения аргумента искомых функций, во втором и третьем столбцах – значения самих функций при соответствующем значении аргумента. При построении графика можно воспользоваться первым столбцом полученной матрицы как аргументом, а вторым и третьим столбцами – как функциями.