Решение дифференциальных уравнений в эксель
Уравнение называется обыкновенным дифференциальным n-го порядка, если F определена и непрерывна в некоторой области и, во всяком случае, зависит от . Его решением является любая функция u(x), которая этому уравнению удовлетворяет при всех x в определённом конечном или бесконечном интервале. Дифференциальное уравнение, разрешенное относительно старшей производной имеет вид
Решением этого уравнения на интервале I=[a,b] называется функция u(x).
Решить дифференциальное уравнение у / =f(x,y) численным методом - это значит для заданной последовательности аргументов х0, х1…, хn и числа у0, не определяя функцию у=F(x), найти такие значения у1, у2,…, уn, что уi=F(xi)(i=1,2,…, n) и F(x0)=y0.
Таким образом, численные методы позволяют вместо нахождения функции y=F(x) (3) получить таблицу значений этой функции для заданной последовательности аргументов. Величина h=xk-xk-1 называется шагом интегрирования.
Метод Эйлера относиться к численным методам, дающим решение в виде таблицы приближенных значений искомой функции у(х). Он является сравнительно грубым и применяется в основном для ориентировочных расчетов. Однако идеи, положенные в основу метода Эйлера, являются исходными для ряда других методов.
Метод Эйлера для обыкновенных дифференциальных уравнений используется для решений многих задач естествознания в качестве математической модели. Например задачи электродинамики системы взаимодействующих тел (в модели материальных точек), задачи химической кинетики, электрических цепей. Ряд важных уравнений в частных производных в случаях, допускающих разделение переменных, приводит к задачам для обыкновенных дифференциальных уравнений – это, как правило, краевые задачи (задачи о собственных колебаниях упругих балок и пластин, определение спектра собственных значений энергии частицы в сферически симметричных полях и многое другое)
Обзор методов решения в Excel
1.1 Метод Рунге-Кутта четвертого порядка для решения уравнения первого порядка
Идея Рунге-Кута состоит в том, чтобы использовать метод неопределённых коэффициентов. Наиболее употребительным методом Рунге-Кутта решения уравнения первого порядка y' = F(x,y) (1) является метод четвертого порядка, в котором вычисления производятся по формуле:
Расчетная таблица, скопированная с рабочего листа MS Excel, показана на рис. 2.1.
Рис. 2.1. Фрагмент рабочего листа с решением примера методом Эйлера
2.1.4.Численное решение методом Рунге-Кутта в Excel
Расчетные формулы для конкретного примера запишем, исходя из общих формул и подставляя в них конкретные функции, определяемые видом уравнения (2.1). Отметим, что исходное уравнение должно быть приведен к виду (1.1).
Задаемся начальным значением y0 и вычисляем
Приведем также пример расчетных формул для уравнения другого вида, а именно уравнения, содержащего x в правой части.
Для него расчетные формулы имеют вид:
Расчетная таблица, скопированная из Excel, представлена на рис. 2.2. Ниже указан порядок ее заполнения
Первая строка заполнена именами переменных
Две ячейки H21, H22 отводятся под значения констант h, c и еще две ячейки G21?G22 под их имена
Первый столбец заполняется значениями x
В ячейку B22 вводим значение y0,
В ячейки С22,D22, E22, F22 вводим соответственно формулы:
В ячейку B23 вводим формулу:
В ячейки B24:B37 копируем формулу из ячейки B23
В ячейки C22:F36 копируем формулы из ячеек C22:F22.
Рис. 2.2. Фрагмент рабочего листа с решением примера по методу Рунге-Кутта
Построим графики численного решения данного уравнения в Excel методами Эйлера и Рунге-Кутта, а также график аналитического решения. Все три графика строим на одной диаграмме. При построении графика с помощью Мастера диаграмм следует:
1. На первом шаге выбрать Тип диаграммы - Точечная
2. На втором шаге последовательно задать ряды данных для каждого графика:
3. На третьем шаге задать параметры Диаграммы: заголовок, названия осей. В случае абстрактного уравнения заголовок диаграммы можно опустить. В легенде следует указать, какому методу данный ряд соответствует.
4. На последнем шаге помещаем диаграмму на имеющемся рабочем листе.
1. Если графики трудно различимы, можно изменить формат рядов данных (на графике), выделяя соответствующий ряд и вызывая контекстное меню правой кнопкой мыши
2. Если шкалы по осям x или y не соответствуют диапазону изменения x, y, нужно изменить их минимальные и максимальные значения по команде формат оси, вызываемой в контекстном меню.
Графики решения данного уравнения различными методами представлены на рис. 2.3.
Рис. 2.3. Графики решений уравнения (2.1) различными способами
2.1.5. Численное решение задачи методом Рунге-Кутта в Турбо Паскале
Расчетные формулы.
Решим уравнение (2.1) методом Рунге-Кутта четвертого порядка, используя формулы (1.8)-(1.9). Т.к. неизвестные значения функции y(x) вычисляются рекуррентно, требуется задать ее начальное значение.
ячейках от C3 до C8 вычисляем значения вспомогательной функции 1/(2 хi + 1), входящей взнаменатели функций p(x), q(x) и f(x). Вводим в C3 формулу "=1/(2*B3+1)" и протягиваем эту формулу до ячейкиC8;
5.
В ячейки D3, E3 и F3 записываем формулы, соответствующие (13), для вычисления значений функций p(x), q(x) и f(x). Запись этих формул при вводе их в ячейки таблицы имеет следующий вид: "=-2*C3", "= ‑12*C3*C3" и "=(3*B3+1)*C3*C3" соответственно. При протягивании этих формул по столбцам D, Е и F до восьмой строки таблица заполняется так, как показано на рис. 2.
6. Далее начинаем заполнение столбцов G, H, I и J значениями коэффициентов Ai, Ci, Bi и Fi в соответствии с форматом системы уравнений (11). В ячейки G3, H3, I3 записываем значения, определяемые форматом первого уравнения системы (11): A0=0, C0=-1, В0=0. В ячейку J3 записываем ссылку на ячейку J1, в которой записано начальное значение F0=Y0: "=J1".
7. В ячейки G8, H8, I8 записываем значения, определяемые конечными условиями A5=0, C5=-1, В5=0. В ячейку J8 записываем ссылку на ячейку L1, в которой записано начальное значение Fk= F5 = Yk :"=L1".
8.
Далее заполняем столбцы G, H, I и J значениями коэффициентов Ai , Ci , Bi и Fi. В ячейку G4 вводим формулу "=1-D4*$H$1/2", соответствующую формуле
для вычисления коэффициента A1. После чего протягиваем эту формулу до ячейки G7.
9.
Аналогично в ячейку Н4 вводим формулу для вычисления коэффициента С1: "2-Е4*$H$1*$H$1", а в ячейку I4 вводим формулу для вычисления коэффициента B1: "1+D4*$H$1/2", реализуя соответствующие формулы
Протягиваем эти формулы до ячеек Н7 иI7.
10.
В ячейках столбца J формируем вектор правых частей системы уравнений (11). В ячейку J4 вводим формулу “=F4*$H$1*$H$1”, соответствующую формуле Fi = fih 2 . Протягиваем эту формулу до ячейки J7. В результате получаем таблицу, показанную на рис. 3. Следует отметить, что в столбцах G, H, I и J этой таблицы записаны элементы матрицы, решаемой системы уравнений (11).
11. Используя вычисленные значения коэффициентов Ai, Ci, Bi и Fi, находим в соответствии с формулами (18) и (19) значения коэффициентов ai и bi. В ячейку K3 запишем формулу для вычисления a0 : "=I3/H3", а в ячейку L3 формулу для вычисления b0.: "=-J3/H3". И далее в ячейки K4 и L4 вводим формулы, соответствующие (19), для вычисления коэффициентов
a1 :"=I4/(H4-K3*G4)", и b1: "=(G4*L3‑J4)/( H4‑K3*G4)".
Протягиваем эти формулы до ячеек K8 и L8 соответственно. Результаты вычисления показаны на рис. 4.
12.
Дальнейшие вычисления выполняются в столбце М по формулам обратного хода (21) и (22). В ячейку М8 введем формулу "=L8", представляющую ссылку на значение bn. В ячейку М7 введем формулу "=L7+K7*M8", соответствующую и протянем ее до ячейки М3. Результаты вычислений показаны на рис. 5.
13. Построим график функции Y(X), используя возможности мастера диаграмм программы MS Excel. Для этого выделим ячейки столбца значений функции Y(X) от ячейки М3 до М8 ивыполним процедуру создания диаграммы, используя средства "мастера диаграмм" программы MS Excel. Окончательный результат показан на рис. 5.
Проверка правильности полученного решения
Для проверки правильности полученного в столбце М решения Y(X) выполняется подстановкой полученного решения в уравнения исходной системы. В ячейки столбца N последовательно вводим формулы, реализующие вычислительный алгоритм, определяемый системой уравнений (11). В N3 введём формулу "=M3". Затемв N4 введём формулу "=G4*M3‑H4*M4+I4*M5", реализующую левую часть уравнения (9). Эту формулу протягиваем до N8, получаясоответственно – "=G5*M4‑H5*M5+I5*M6" вN5, "=G6*M5-H6*M6+I6*M7" в N6, "=G7*M6‑H7*M7+I7*M8" вN7 и "=G8*M7-H8*M8+I8*M9" вN8. Полученные в столбце N значения совпадают со значениями в столбце J, что позволяет судить о правильности полученного решения.
Варианты заданий для выполнения самостоятельной работы
Методом конечных разностей найти решение краевой задачи на сетке из 6 узлов.
Вариант 9
Вариант 10
Вариант 11
Вариант 12
Вариант 13
Вариант 14
Вариант 15
Литература
1. Пискунов Н.С. Дифференциальное и интегральное исчисления. М.:Физматгиз, 1963. – 856 с.
2. Карпов В.В., Коробейников А.В. Математические модели задач строительного профиля и численные методы их исследования: Учеб. пособие. – СПб., СПбГАСУ, 1996. – 134 с.
3. Бахвалов Н.С. Численные методы. – М.:Наука, 1973. – 632 с.
4. Вагер Б.Г. Численные методы решения дифференциальных уравнений: Учеб. пособие – СПбГАСУ. – СПб. 2003. 114 с.
5. Любимов Е.Б. и др. Решение систем линейных алгебраических уравнений средствами программы Microsoft Excel: Метод. указ. – СПб., СПбГАСУ, 2005. – 22 c.
Основные понятия, используемые в постановках краевых задач. 3
Применение метода прогонки для решения систем линейных алгебраических уравнений с трёхдиагональными ленточными матрицами. 6
Реализация метода прогонки в среде программы MS Excel 6
Постановка задачи. 6
Проверка правильности полученного решения. 12
Варианты заданий для выполнения самостоятельной работы.. 13
[1] ) Формулы вводятся в ячейки таблиц, начиная с символа “=” (равно). Двойные кавычки использованы в тексте для выделения формулы. Вводить их в ячейки таблицы не нужно.
[2] ) Терминология и сокращения, используемые в тексте методических указаний, приведены в начальном разделе методических указаний к первой лабораторной
Данная система может решаться как система двух уравнений первого порядка (1.10)-(1.11), где f1 (x,y, z)=z , f2 (x,y,z)=f(x, y, z) .
2. Примеры постановок и численного решения задач.
2.1. Модель размножения бактерий ( решение задачи Коши для уравнения первого порядка)
2.1.1. Постановка задачи.
Рассмотрим задачу размножения бактерий при условии, что скорость прироста числа бактерий пропорциональна их наличному количеству с коэффициентом пропорциональности с . Задаваясь начальным количеством бактерий y(0)=y0 в начальный момент времени t=0 и некоторым промежутком времени T, определить закон изменения числа бактерий на промежутке времени [0,T] . Задачу решить численно методами Эйлера и Рунге-Кутта с использованием табличного процессора MS Excel и системы программирования Турбо Паскаль 7.0.
При выполнении расчетов положить y(0)=y0, c=1, h=0.1, T=1.5.
2.1.2. Математическая модель задачи
Пусть y=y(t) - количество бактерий в произвольный момент времени. Скорость роста бактерий представляет собой производную функции y(t). По условию выполнено соотношение:
Полагая y(0)=y0 , получим задачу Коши для дифференциального уравнения первого порядка.
2.1.3. Численное решение методом Эйлера
Расчетные формулы.
Данное уравнение имеет аналитическое решение
что предоставляет возможность сравнить с ним численное решение, полученное тем или иным методом.
При выполнении расчетов в Excel целесообразно записать формулу метода Эйлера (1.4) с указанием конкретной правой части уравнения.
Для проведения расчетов будем вычислять вначале добавку к текущему значению функции для вычисления следующего
Выполнение расчетов в Excel
Рассмотрим расчетную таблицу в Excel, содержащую три столбца для значений . Дадим им заголовки x, y, , расположив их в ячейках A2:C2. Для постоянных величин h, c отведем отдельные ячейки E2 и E3. Их заголовки помещены в ячейки D2 и D3.
Расчеты в таблице Excel выполняются по следующему алгоритму
1. Вычисление первого столбца:
первые два значения x = x0 и x1 = x0+h вводятся в ячейки A3 и A4, затем, выделив две эти ячейки , заполняем столбец значений x до достижения конечного значения x=T.
2. Затем заполним первую строку расчетной таблицы: в столбце для y введем y0 в ячейку B3, в столбце для введем в ячейке C3 формулу:
=$E$2*$E$3*B3 (вычисляется приращение функции y для текущего значения x в соответствии с по формуле 2.4)
3. Вычисление второго столбца:
Вводим формулу =В3+С3 в ячейку B4 и копируем ее в ячейки B5:B18 (вычисляется новое значение функции y при изменении x на один шаг с помощью линейного приращения по формуле 2.5).
Читайте также: