Численные методы в среде символьной математики. Программирование и решение задач: Учебно-методическое пособие, страница 37

Методом Рунге-Кутта четвертого порядка проинтегрировать систему линейных дифференциальных уравнений, полученную в задаче № 16, которые описывают остывание неравномерно нагретого по длине теплоизолированного стержня. Масштабный коэффициент В взять равным единице, шаг по безразмерной пространственной переменной  установить в соответствии с выбранной многоточечной аппроксимацией производных. Кривые распределения температур в заданных точках стержня изобразить на графике  в моменты времени  секунд.

Подготовка параметров процедуры Рунге-Кутта

Построение процедуры для численного интегрирования одного дифференциального уравнения первого порядка методом Рунге-Кутта было продемонстрировано в задаче №21. В данной задаче требуется проинтегрировать систему уравнений первого порядка. Чтобы ускорить работу по получению таблицы решений системы, применим имеющуюся в утилите ODE_APPR.MTH  процедуру Рунге-Кутта.

Предварительно вычислим коэффициент k в системе уравнений и оценим величину шага по времени ht:

Опишим каждое уравнение системы, вспомнив, что на границах ¾0=¾7=0.

u1:=k·(-20·¾1+6·¾2+4·¾3-¾4)

u2:=k·(16·¾1-30·¾2+16·¾3-¾4)

u3:=k·(-¾1+16·¾2-30·¾3+16·¾4-¾5)

u4:=k·(-¾2+16·¾3-30·¾4+16·¾5-¾6)

u5:=k·(-¾3+16·¾4-30·¾5+16·¾6)

u6:=k·(-¾3+4·¾4+6·¾5-20·¾6)

Вектор из уравнений, представляющий систему:

u:=[u1,u2,u3,u4,u5,u6]

Сформируем вектор переменных, входящих в систему, и вектор их начальных значений:

v:=[t,¾1,¾2,¾3,¾4,¾5,¾6]

v0:=[0,-0.06,0.664,0.252,-0.372,-1.424,-1.94]

Так как по условию задачи значения функции в точках необходимо выводить в моменты времени [0. 0.1, 0.2, 0.6, 1, 2, 5] секунд, то шаг по времени необходимо выбрать кратным этим значениям, но с учетом неравенства, приведенного выше. Ближайшим является ht=0.02 с.

С учетом оговоренного процедуру Рунге-Кутта запишим так:

RK(u,v,v0,ht,n)

Здесь неопределенным пока является число итераций n. Чтобы при шаге по времени 0.02 с получить отсчет в момент времени, равный 5 с, необходимо выполнить 250 рекуррентных подстановок. Такую таблицу выводить на экран крайне неразумно, поэтому из нее, как из вектора с векторными компонентами, сделаем выборку по индексу, представленному списком значений : [1. 6, 11, 31, 51, 101, 251]

NotationDigits:=4

VECTOR((RK(u,v,v0,0.02,250)) ,i,[1,6,11,31,51,101,251])

                            i                          

„  0      -0.06        0.664        0.252       -0.372       -1.424        -1.94    †

¦                                                                                   ¦

¦ 0.1    0.02980     -0.006536     -0.1812      -0.4400      -0.5917      -0.4419   ¦

¦                                                                                   ¦

¦ 0.2   -0.05162      -0.1237      -0.2148      -0.2869      -0.2856      -0.1807   ¦

¦                                                                                   ¦

¦ 0.6   -0.03062     -0.05546     -0.06957     -0.06998     -0.05640     -0.03136   ¦

¦                                                                                   ¦

¦  1    -0.008523    -0.01538     -0.01919     -0.01919     -0.01538     -0.008528  ¦

¦                                                                                   ¦

¦  2   -0.0003382   -0.0006103   -0.0007614   -0.0007614   -0.0006103   -0.0003382  ¦

¦                                                                                   ¦

¦               -8           -8           -8           -8           -8           -8 ¦

…  5   -2.112·10    -3.811·10    -4.755·10    -4.755·10    -3.811·10    -2.112·10   ‡

Как видим процесс охлаждения стержня практически завершился в момент времени, равный 1 секунде. Поэтому графики распределения температур по длине стержня изобразим лишь для первых пяти моментов времени: [0. 0.1, 0.2, 0.6, 1].