Численное решение ОДУ в Mathcad: лабораторная работа

ЛАБОРАТОРНАЯ РАБОТА № 4
Тема работы: Численное решение обыкновенных дифференциальных уравнений. Оценка
погрешности решения. Визуально-ориентированное и модульное программирование в
среде Mathcad.
Цель работы: Усвоить алгоритмы классических численных методов интегрирования
обыкновенных дифференциальных уравнений (Рунге-Кутта, Эйлера). Освоить технику
создания законченных программных модулей (блоков) в среде Mathcad.
Задание: Задана система обыкновенных дифференциальных уравнений (варианты
заданий совпадают с лабораторной работой №3). Численно решите задачу Коши для
данной системы уравнений:
- используя встроенную функцию rkfixed системы Mathcad;
используя
средства
визуально-ориентированного
программирования
(ранжированные переменные) и методы Рунге-Кутта и Эйлера;
- используя средства модульного программирования и метод Рунге-Кутта.
Сравните численные результаты между собой. Оцените погрешность численного
интегрирования по правилу Рунге.
Теоретическая часть
Рассмотрим задачу Коши для одного дифференциального уравнения первого порядка
du
u t t  u 0 .
 f t ,u t  ,
0
dt
Для численного решения этого уравнения зададим сетку t k  t 0  k, k  0,1, , N  с
постоянным шагом  . На этой сетке определим сеточные функции yk  u tk  ,
f k  f tk ,uk  .
Метод Рунге-Кутта позволяет строить схемы различного порядка точности.
Наиболее известной и широко используемой на практике является схема Рунге-Кутта
четвертого порядка точности:
k1  f tn , yn 

 

k 2  f  t n  , yn  k1 
2
2 


 

k3  f  t n  , yn  k 2 
2
2 

k4  f tn  , yn  k3 

yn 1  yn  k1  2k2  2k3  k4  .
6
Численный метод Эйлера
y n1  y n  f t n , y n 
является частным вариантом метода Рунге-Кутта первого порядка точности.
Метод Рунге-Кутта легко обобщается на системы уравнений путем формальной
замены скалярных величин u , f t ,u  на векторы u ,f t ,u  . Для системы уравнений
du
T
T
 f t ,ut , u  u1 ,u 2 , ,u m  , f   f1 , f 2 , , f m 
dt
расчетные формулы имеют вид:
k 1  f t n , y n 

 

k 2  f  t n  ,y n  k 1 
2
2 


 

k 3  f  t n  ,y n  k 2 
2
2 

k 4  f t n  , y n  k 3 

y n1  y n  k 1  2k 2  2k 3  k 4 
6
В частности, для системы двух дифференциальных уравнений
 du
 dt  f t ,u t , vt 
 dv
  g t ,u t , vt 
 dt
развернутая запись формул схемы Рунге-Кутта имеет вид:
k1  f t n , y n , z n  ,
q1  g t n , y n , z n  ,


 


 


k 2  f  t n  , y n  k1 , z n  q1  ,
q 2  g  t n  , y n  k1 , z n  q1  ,
2
2
2 
2
2
2 




 


 


k3  f  t n  , yn  k 2 , z n  q2  ,
q3  g  t n  , y n  k 2 , z n  q 2  ,
2
2
2 
2
2
2 


k 4  f t n  , y n  k 3 , z n  q3  ,
q4  g t n  , y n  k 3 , z n  q3  ,


y n 1  y n  k1  2k 2  2k 3  k 4  ,
z n1  z n  q1  2q 2  2q3  q 4 
6
6
Здесь через y , z обозначены приближенные сеточные функции, соответствующие
функциям ut ,vt  соответственно.
Шаг интегрирования в методах Эйлера и Рунге-Кутта следует выбирать достаточно
малым, чтобы обеспечить устойчивость и требуемую точность расчета. Как практически
получить нужную точность? Априорные оценки точности для выбора шага  на деле не
используются по следующим причинам. Во-первых, они чрезвычайно громоздки и
включают производные, которые до начала расчета не известны, во- вторых, эти оценки
являются мажорантными и могут во много раз превосходить фактическую ошибку
расчета. Поэтому основным практическим приемом является апостериорная оценка
точности. Для ее получения расчет проводят на двух сетках с шагом  и k , k  1 , затем
применяют правило Рунге. Если численные решения на двух сетках обозначить через
yt ;  и yt ; k , то погрешность решения на сетке с меньшим шагом оценивается по
формуле
y t ;   y t ; k
t ;  
,
k p 1
где p - порядок точности схемы.
Порядок выполнения задания
1. Определите функции, входящие в правую часть системы уравнений.
2. Задайте исходные данные для решения задачи Коши, и вычислите шаг
интегрирования.
3. Используя ранжированные переменные, составьте программу численного решения
задачи Коши методом Рунге-Кутта и методом Эйлера.
4. Вычислите решение методами Рунге-Кутта и Эйлера.
5. Сравните полученные решения, изобразив их на одном графике.
6. Решите задачу Коши с использованием встроенной функции rkfixed.
7. Убедитесь, что это решение совпадает с полученным ранее (в п.4) решением.
8. Увеличив в два раза число узлов интегрирования, и заново решив задачу Коши,
оцените погрешность численного решения.
9. Используя средства модульного программирования системы Mathcad, составьте
законченный модуль численного решения задачи Коши для системы двух
дифференциальных уравнений методом Рунге-Кутта с постоянным шагом  .
10. Убедитесь в правильной работе данного модуля, сравнив полученное с его
помощью решение с найденным ранее (п.4, п.6).
11. Сформулируйте выводы по проделанной работе.
12. Сохраните рабочий документ.
В приведенном ниже примере предполагается, что система дифференциальных уравнений
имеет вид:

 du1
2
2
 dt  u 2  u1 4  u1  u 2
 du
 2  u1  u 2 4  u12  u 22
 dt


