CC BY
5
Discrete & Continuous Models & Applied Computational Science 2023, 31 (2) 164-173
ISSN 2658-7149 (online), 2658-4670 (print) http://journals-rudn-ru/miph
Research article
UDC 519.872:519.217 PACS 07.05.Tp, 02.60.Pn, 02.70.Bf DOI: 10.22363/2658-4670-2023-31-2-164-173 EDN: XDVQBB
Implementation of the Adams method for solving ordinary differential equations in the Sage computer
algebra system
Mikhail D. Malykh1'2, Polina S. Chusovitina1
1 RUDN University, 6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation 2 Meshcheryakov Laboratory of Information Technologies, Joint Institute for Nuclear Research, 6, Joliot-Curie St., Dubna, Moscow Region, 141980, Russian Federation
(received: April 25, 2023; revised: May 7, 2023; accepted: June 26, 2023)
Abstract. This work is devoted to the implementation and testing of the Adams method for solving ordinary differential equations in the Sage computer algebra system. The Sage computer algebra system has, to some extent, trivial means for numerical integration of ordinary differential equations, but at the same time, it is worth noting that this environment is convenient and practical for conducting computer experiments related to symbolic numerical calculations in it. The article presents the FDM package developed on the basis of the RUDN, which contains the developments of recent years, performed by M.D. Malykh and his students, for numerical integration of differential equations. In this package, attention is paid to the visualization of the calculation results, including the construction of various kinds of auxiliary diagrams, such as Richardson diagrams, as well as graphs of dependence, for example, the value of a function or step from a moment in time. The implementation of the Adams method will be considered from this package. In this article, this implementation of the Adams method will be tested on various examples of input data, and the method will also be compared with the Jacobi system. Exact and approximate values will be found and compared, and an estimate for the error will be obtained.
Key words and phrases: differential equations, Adams method, Sage, FDM package, Cauchy theorem, Taylor series, Richardson diagram
1. Introduction
To describe models in a variety of subject areas from mechanics to economics, ordinary differential equations are used [1]. These equations admit solutions in elementary functions only in some very special cases, therefore they are usually
© Malykh M.D., Chusovitina P.S., 2023
This work is licensed under a Creative Commons Attribution 4.0 International License https://creativecommons.org/licenses/by-nc/4-0/legalcode
investigated numerically. The finite difference method was proposed by Euler, the Runge-Kutta method of the 4th order is the most popular numerical method for solving initial problems for ordinary differential equations [2].
Old authors, including J. Scarborough [3, ch. XIII], mention numerical methods alternative to the Runge-Kutta method. The method that J. Scarborough has associated with the name of the English theoretical astronomer J. K. Adams, was forgotten for a long time, because it was very inconvenient to implement on a computer: before its use, a number of preparatory calculations had to be carried out on paper. However, with the development of computer computing, it became possible to perform these actions on a computer, which pushes us to study the possibility of implementing the Adams method in modern computer algebra systems.
Currently, RUDN University is developing an addition to Sage — the FDM package, which contains the achievements of recent years, made by M. D. Malykh and his students. The goal of the project is to create a convenient environment for numerical experiments with ODES in the Sage computer algebra system. This project is available to everyone on https://github. com/malykhmd/fdm. The general principles of the package are described in [4].
The purpose of this work is to test the implementation of the Adams method in FDM.
2. The Adams method and its implementation in FDM
Consider the initial problem
dit=fM> x(0) = xo. (1)
Its solution exists by virtue of the Cauchy theorem. Decompose its solution into a Taylor series at dt = 0:
x(t + dt) = x(t) + x(t)dt + 1x(t)dt2 + .... (2)
Scarborough was forced to search numerically for the coefficients of the Taylor series, but now we can find them analytically, using the formulas specified by Cauchy himself. If f(x, t) is known, then
x(t) = f(x(t),t), x = Df, x = D2 f,
where
Dg = f— + — dx dt'
To calculate the coefficients of the Taylor series using these formulas, it is required to differentiate the symbolic expression f many times, which is naturally performed in Sage.
The Adams method in our interpretation is as follows. First, according to the given symbolic expression f, the Taylor polynomial is compiled
x + f(x,t)-dt + --- + 1Dr-1 (f)-dtr (3)
T.
up to r-th order members. Its coefficients are calculated explicitly in symbolic form.
Then the segment 0 < t < Ton which the initial problem is considered is divided into segments of width dt. At time t = 0, the solution is given to us. To find a solution at time t = dt, we substitute t = 0, x = xn into the Taylor polynomial (3) and thus obtain an approximate value of xx for x(dt). To find a solution at time t = 2 • dt, we substitute t = dt, x = x^, etc. into the Taylor polynomial.
FDM implements two versions of the Adams method: with a constant step, and with a step that becomes smaller the larger the first discarded term
VTTyD'(f)'dt'+1
in the Taylor formula (3).
3. Numerical experiments
We can test the implementation of the Adams method in FDM with a few examples. We will evaluate the error using the Richardson method [2].
Example 1. Consider the initial problem
= t2 + x, x(0) = 0
dt
on the segment 0 < t < 1, its exact solution is known: var("x,t")
prl=Initial_problem(x,t~2+x,0,1)
P=adams_adaptive(prl, h=0.1, field=RealField(500))
We get the value of x at the point t = 0.8: P.value(x,0.8)
0.21079360329767660164890230589662678539752960205078125 In figure 1 we see that the exact solution coincides with the approximate one. In figure 2 we see that the step depends on the moment of time t, in this case it decreases monotonically.
0.145 0.140 0.135 0.130 0.125 0.120 0.115 0.110 0-105
0.2 0.4 0-6 0.8
Figure 1. Dependence of x on t
Figure 2. Dependence of the step on the moment of time t
Figure 3 shows the log-log plot for the dependence of the error E on the value of the variable x at time t = 0.9 on the parameter h characterizing the step of the Adams method. Such a figure is known as Richardson diagram [2]. Theoretically, the error E = chr + 0(hr+1).
In our experiments we use 3 terms in Taylor series, thus r must be equal 3, so the error must be proportional to h3. In our example (figure 3) there is a direct proportionality on the logarithmic scale and the slope is equal to 2.99, which is close to 3, so the error is proportional to h3. To build a Richardson diagram [2] we used the standard tools of FDM [4]:
L= [adams_adaptive(prl, h=l/10/2~n, field=RealField(500))
for n in range(lO)] richardson_plot(L,x,0.9)
№)l
Figure 3. Richardson diagram
Each element of the list is an approximate solution with respect to the step h = 0.1 • (l/2n), n = 0,1,...,9.
Using standard tools of FDM, we can calculate the values of x and the error and make sure that it is small. For example, at time t = 0.8 and h = 1/20:
richardson(L[0], L[l], x, 0.8)
[0.21104359494629298943, -0.000035713092659491899]
If we take the exact solution and subtract the approximate one from it, we get an error value equal to —0.0000382:
x_exact=-t~2 - 2*t + 2*e~t - 2 L[l].value(x,0.8) - RR (x_exact.subs(t=0.8)) -0.0000382620386425447 It should be noted that the values are almost the same (—0.0000357 and —0.0000382), which means that Richardson's method gives an almost correct estimate for the error.
Example 2. Now we consider the system of two ODEs
/y> - rr* rp - -rr*
•Aj — 2 — X
with initial conditions 1^(0) = 0, ;r2(0) = 1 at the interval 0 < t < 10. Let's take the step h = 1/4 and solve the initial problem in FDM:
var("xl,x2,t")
pr2=Initial_problem([xl,x2] ,[x2,-xl] , [0,1] ,10) P=adams_adaptive(pr2,h=l/4)
In figure 4 we can see that our points lie on the sine wave. To plot graphs we use the line:
P.plot(t,xl,color='black') + plot(sin,(0,10),color='grey')
Xl
Figure 4. Dependence of x on t
To plot a Richardson diagram (figure 5) we use the line: L=[adams_adaptive(pr2,h=l/4/2"n) for n in range(10)]
№i)l
Figure 5. Richardson diagram
In this Richardson diagram we can see that not all points lie on a straight line, this is explained by the deviation of the slope value, which is not 3, but 3.09. But the error is also proportional to h3. Calculate the values of x and the error:
richardson (L[-3],L[-2], x 1,9)
[0.412118484375421, -8.69428719494993e-10]
We can compare the received error value by the Richardson method with the real error:
L[-2].value(xl,9)-sin(9).n() -8.66335503335591e-10
Note again that the values are almost identical.
Example 3. Now consider the Jacobi oscillator
p = qr, q = —pr, r = —k2pq
with initial conditions p{0) = 0, q{0) = 1, r(0) = 1 at the interval 0 < t. < 100. This example is interesting as a nonlinear oscillator. We described the problem by lines:
var("p,q,r,t") k=0.99
pr3=Initial_problem([p,q,r], [q*r,-p*r,-k~2*p*q], [0,1,1], 100) P=adams_adaptive(pr3,h=l/4)
and plot the graphs to compare approximate and exact solutions in Jacobi elliptic functions (figure 6):
P.plot(t,p,color='grey')+plot(jacobi('sn', t, k~2),(t,0,100), color='black')
P
1.0-
A
A
1
r\
0.5-
2D
100
-0.5-
-1.0-
X
\J
\J
V.
Figure 6. Comparison of values obtained by the Adams method and the Jacobi method
It is noticeable that the values diverge over a long-time interval, so that this does not happen, it is necessary to take a smaller step. In figure 7 we use a small-time interval, and the solutions coincide.
Figure 7. Comparison of exact and approximate solutions
-0.5
-1.0
We use adaptive Adams method what is important for nonlinear problems. In figure 8, we see the dependence of the step on the moment of time, while in places where the values change smoothly, the step is larger, and where the function changes quickly, the step is smaller, which indicates that the chosen method of adaptation according to the last term in the Taylor series is correct.
dt
0.8
0.7
0.6
0.5
0.4
0.3
10
Figure 8. Dependence of the step on the moment of time t
At Richardson diagram (figure 9), we can see that the points lie on a straight line with a slope of 2.99. The error is also proportional to h3. Also we calculate the value of x, the error by the Richardson method and the real error, and again we get identical error values: richardson(L[3],L[4],p,9)
[-0.983742403638134, -2.91251782518526e-6] richardson(L[3],L[4],p,9)[0]-jacobi('sn',9,k~2) -2.91556109444091e-6
№)l
Figure 9. Richardson diagram
4. Conclusion
In this article, the Adams method for solving ordinary differential equations was tested on three examples. Also, shortcomings were identified in the fdm.sage package, after their correction, the execution of the examples became correct.
Thus, the computer experiments carried out confirm that the correct approximation method is indicated in the implementation of the Adams method in FDM, and this implementation itself allows to search for solutions with an accuracy close to the rounding error. Also, the achievement of the method is the natural adaptation of the step to changes in the function, where the function changes quickly, the step decreases. In this case, the method is symbolic-numerical, since the Taylor series is calculated symbolically once, and then used as a formula in which specific numeric values are substituted. At the beginning of the 20th century, this was considered the main drawback of the method proposed by Adams, since there were problems with symbolic computation - many iterations of differentiation led to large and complex expressions. And now, on the contrary, we can consider the advantage of the method that it can combine the powers of both symbolic and numerical methods, but, of course, the implementation of this method can still be refined and optimized.
Acknowledgments
This work is supported by the Russian Science Foundation (grant no. 2011-20257).
References
[1] H. Gould and J. Tobochnik, An introduction to computer simulation methods. Applications to physical systems. Addison-Wesley Publishing Company, 1988.
[2] A. Baddour and M. D. Malykh, "Richardson-Kalitkin method in abstract description," Discrete and Continuous Models and Applied Computational Science, vol. 29, no. 3, pp. 271-284, 2021. DOI: 10.22363/2658-46702021-29-3-271-284.
[3] J. B. Scarborough, Numerical methods of mathematical analysis. Oxford book company, 1930.
[4] L. Gonzalez and M. D. Malykh, "On a new package for numerical solution of ordinary differential equations in Sage [O novom pakete dlya chislennogo resheniya obyknovennykh differentsial'nykh uravneniy v Sage],"in Proceedings of ITTMM'22, Moscow, Russia, in Russian, 2022, pp. 360-364.
For citation:
M. D. Malykh, P. S. Chusovitina, Implementation of the Adams method for solving ordinary differential equations in the Sage computer algebra system, Discrete and Continuous Models and Applied Computational Science 31 (2) (2023) 164-173. DOI: 10.22363/2658-4670-2023-31-2-164-173.
Information about the authors:
Malykh, Mikhail D. — Doctor of Physical and Mathematical Sciences, Assistant Professor of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia named after Patrice Lumumba (RUDN University) (e-mail: malykh-md@rudn.ru, phone: +7(495)9550927, ORCID: https://orcid.org/0000-0001-6541-6603, ResearcherlD: P-8123-2016, Scopus Author ID: 6602318510)
Chusovitina, Polina S. — Student of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia named after Patrice Lumumba (RUDN University) (e-mail: 1032192941@rudn.ru, ORCID: https://orcid.org/0009-0006-4191-2454)
УДК 519.872:519.217 PACS 07.05.Tp, 02.60.Pn, 02.70.Bf DOI: 10.22363/2658-4670-2023-31-2-164-173 EDN: XDVQBB
Реализация метода Адамса для решения обычных дифференциальных уравнений в системе компьютерной алгебры Sage
М. Д. Малых1'2, П. С. Чусовитина1
1 Российский университет дружбы народов, ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия 2 Лаборатория информационных технологий им. М. Г. Мещерякова,
Объединённый институт ядерных исследований, ул. Жолио-Кюри, д. 6, Дубна, Московская область, 141980, Россия
Аннотация. Работа посвящена реализации и тестированию метода Адамса для решения обыкновенных дифференциальных уравнений в системе компьютерной алгебры Sage. Система компьютерной алгебры Sage обладает в какой-то степени тривиальными средствами для численного интегрирования обыкновенных дифференциальных уравнений, но при этом, стоит заметить, что данная среда удобна и практична для проведения в ней компьютерных экспериментов, связанных с символьно-численными вычислениями. В работе представлен пакет FDM, разработанный на базе РУДН, в котором собраны наработки последних лет, выполненных М. Д. Малых и его учениками, для численного интегрирования дифференциальных уравнений. В данном пакете уделено внимание визуализации результатов вычисления, в том числе построению разного рода вспомогательных диаграмм, например диаграмм Ричардсона, а также графиков зависимости, например значения функции или шага от момента времени. В статье рассмотрена реализация метода Адамса, проведено её тестирование на различных примерах входных данных, а также выполнено сравнение метода с системой Якоби. Найдены и точные, и приближённые значения, проведено их сравнение, получена оценка для ошибки.
Ключевые слова: дифференциальные уравнения, метод Адамса, Sage, пакет FDM, теорема Коши, ряд Тейлора, диаграмма Ричардсона
