# Generalized Runge-Kutta Method with respect to the Non-Newtonian Calculus

Authors: Uğur Kadak; Muharrem Özlük
Journal: Abstract and Applied Analysis
 (2015)

---

## Abstract

Theory and applications of non-Newtonian calculus have been evolving rapidly over the recent years. As numerical methods have a wide range of applications in science and engineering, the idea of the design of such numerical methods based on non-Newtonian calculus is self-evident. In this paper, the well-known Runge-Kutta method for ordinary differential equations is developed in the frameworks of non-Newtonian calculus given in generalized form and then tested for different generating functions. The efficiency of the proposed non-Newtonian Euler and Runge-Kutta methods is exposed by examples, and the results are compared with the exact solutions.


---

## Body

## 1. Introduction

Differential and integral calculus, the most applicable mathematical theory, was created independently by Newton and Leibnitz in the second half of the 17th century. Later Euler redirected calculus by giving a central place to the concept of function and thus founded analysis. Two operations, differentiation and integration, are basic in calculus and analysis. In fact, they are the infinitesimal versions of the subtraction and addition operations on numbers, respectively. The invention of multiplicative calculus can be dated back to 1938, when Volterra and Hostinsky [1] proposed the so-called Volterra calculus that was later identified as a special case of multiplicative calculus. In the period from 1967 till 1972, Grossman and Katz [2] introduced the non-Newtonian calculus consisting of the branches of geometric, bigeometric, quadratic, and biquadratic calculus and so forth. Also Grossman extended this notion to the other fields in [3, 4]. All these calculi can be described simultaneously within the framework of a general theory. We prefer to use the namenon-Newtonian to indicate any of calculi other than the classical calculus.Many authors have extensively developed the notion of multiplicative calculus [5–8]. The complete mathematical description of multiplicative calculus was given by Bashirov et al. [9, 10]. Based on this work various multiplicative numerical approximation techniques were proposed and discussed. Multiplicative calculus has also found its way into biomedical image analysis or modelling with differential equations [11]. Furthermore even the Runge-Kutta method was developed in the framework of bigeometric or Volterra calculus for applications in dynamic systems by Aniszewska [12]. In some cases, for example, for wage-rate (in dollars, euro, etc.) related problems, the use of bigeometric calculus is advocated instead of a traditional Newtonian one. Multiplicative calculus has a relatively restrictive area of applications compared to the calculus of Newton and Leibnitz. Indeed, it covers only positive functions. In order to circumvent the restriction to positive valued functions of real variable, the geometric multiplicative calculus was extended to complex multiplicative calculus. After the heuristic extension suggested by Uzer [13] a comprehensive mathematical description of the multiplicative complex analysis was presented by Bashirov and Riza [14]. Misirli and Gurefe have introduced multiplicative Adams-Bashforth-Moulton methods for differential equations by numerical examples in [15]. Also some authors have also worked on the classical sequence spaces and related topics by using non-Newtonian calculus [16–18]. Further Kadak [19] and Kadak and Efe [20, 21] have determined Kothe-Toeplitz duals and matrix transformations between certain sequence spaces over the non-Newtonian complex field.It is shown that depending on the problem the non-Newtonian Runge-Kutta method gives satisfactory results according to classical Runge-Kutta method. One of the advantages of this method is that the approximations vary based on the choice of generating functions. Further, the numerical results can be computed by relatively small number of steps for small or large intervals.The rest of this paper is organized as follows. In Section2, we give some required definitions and consequences related to the non-Newtonian calculus, “q” and “exp” generators in particular. In Section 3, first, the definitions of ∗-limit and ∗-continuity are given. Then we construct the notion of ∗-derivative with respect to two generating functions and formulate some useful results of differential calculus in terms of ∗-derivative which will be used in the proof of our main results. Also higher-order ∗-derivatives of a function are given. In Section 4, first we give Lipschitz condition for the existence and uniqueness for an initial value problem. Non-Newtonian Taylor expansions of functions of several variables are presented according to classical Taylor expansion. Secondly we derive the generalizations of non-Newtonian Euler and Runge-Kutta methods explicitly. Finally two illustrated examples for different cases are analysed and comparison of the numerical results with respect to exact values is given in tables.


## 2. Preliminaries, Background, and Notation

Agenerator is a one-to-one function whose domain is R and whose range is a subset of A
⊂
R. Each generator generates exactly one arithmetic, and conversely each arithmetic is generated by exactly one generator. If I
(
x
)
=
x for all x
∈
R, then I is calledidentity function whose inverse is itself. In the special cases α
=
I and α
=
e
x
p, α generates the classical and geometric arithmetics, respectively. By α-arithmetic, we mean the arithmetic whose domain is R and whose operations are defined as follows: for x
,
y
∈
A
⊂
R and any generator α, (1)
α
-addition
x
+
˙
y
=
α
α
-
1
x
+
α
-
1
y
,
α
-subtraction
x
-
˙
y
=
α
α
-
1
x
-
α
-
1
y
,
α
-multiplication
x
×
˙
y
=
α
α
-
1
x
×
α
-
1
y
,
α
-division
x
/
˙
y
=
α
α
-
1
x
÷
α
-
1
y
,
α
-order
x
<
˙
y
⟺
α
-
1
x
<
α
-
1
y
.As an example, if we choose the function e
x
p from R to the set R
+, (2)
α
:
R
⟶
R
+
x
⟼
y
=
α
x
=
e
x
,and α-arithmetic turns out to be geometric arithmetic: (3)
α
-addition
x
+
˙
y
=
e
l
n
x
+
l
n
y
=
x
·
y
,
α
-subtraction
x
-
˙
y
=
e
l
n
x
-
l
n
y
=
x
÷
y
,
α
-multiplication
x
×
˙
y
=
e
l
n
x
l
n
y
=
x
l
n
y
=
y
l
n
x
,
α
-division
x
/
˙
y
=
e
l
n
x
/
l
n
y
=
x
1
/
l
n
y
.Now we give the infinitely manyq-arithmetics, of which the quadratic and harmonic arithmetics are special cases. The function q
(
x
) and its inverse q
-
1
(
x
) for x
∈
R are defined as follows: (4)
q
x
=
x
1
/
p
,
x
>
0
0
,
x
=
0
-
-
x
1
/
p
,
x
<
0
,
q
-
1
x
=
x
p
,
x
>
0
0
,
x
=
0
-
-
x
p
,
x
<
0
.If p
=
1, then the non-Newtonian calculus is reduced to the classical calculus. Also if p
=
2, then the function q
(
x
) generates the quadratic arithmetic: (5)
q
-addition
x
+
˙
y
=
x
2
+
y
2
1
/
2
=
x
2
+
y
2
,
q
-subtraction
x
-
˙
y
=
x
2
-
y
2
1
/
2
=
x
2
-
y
2
,
q
-multiplication
x
×
˙
y
=
x
2
y
2
1
/
2
=
x
y
,
q
-division
x
/
˙
y
=
x
2
y
2
1
/
2
=
x
y
,
q
-order
x
<
˙
y
⟺
x
1
/
2
<
y
1
/
2
=
x
<
y
.Arithmetic is any system that satisfies the whole of the ordered field axioms whose domain is a subset of R. There are infinitely many types of arithmetic, all of which are isomorphic, that is, structurally equivalent. Nevertheless, the fact that two systems are isomorphic does not exclude their separate usage.


## 3. Definition of the Non-Newtonian Derivative (∗-Derivative)

In the present section we give a new type of calculus denoted by∗-calculus which represents general structure of non-Newtonian calculus. Since all arithmetics are isomorphic, one can easily obtain all arithmetics by using a unique function from a
l
p
h
a to the b
e
t
a arithmetic.Letα and β be two arbitrarily selected generators and let α-arithmetic and β-arithmetic also be the ordered pair of arithmetics. The set (
B
,
+
¨
,
-
¨
,
×
¨
,
/
¨
) is a complete ordered field and β-generator generates β-arithmetic [2]. Definitions given for α-arithmetic are also valid for β-arithmetic.In∗-calculus, α-arithmetic is used for arguments and β-arithmetic is used for values; in particular, changes in arguments and values are measured by α-differences and β-differences, respectively. The operators of this type of calculus are applied only to functions with arguments in A and values in B. The ∗-limit of a function with two generators α and β is defined by (6)
l
i
m
x
→
a
∗
f
x
=
b
⟺
∀
ε
>
¨
0
,
∃
δ
>
˙
0
∋
f
x
-
¨
b
β
≤
¨
ε
,
∀
x
∈
A
,
x
-
˙
a
α
<
˙
δ
.A function f is ∗-continuous at a point a in A if and only if a is an argument of f and ∗
l
i
m
x
→
a
f
(
x
)
=
f
(
a
). When α and β are the identity function I, the concepts of ∗-limit and ∗-continuity are identical with those of classical limit and classical continuity (see [18]).The isomorphism fromα-arithmetic to β-arithmetic is the unique function ι (iota) that possesses the following three properties.(i)
ι is one to one.
(ii)
ι is from A onto B.
(iii)
For any numbersu and v in A, (7)
ι
u
+
˙
v
=
ι
u
+
¨
ι
v
,
ι
u
-
˙
v
=
ι
u
-
¨
ι
v
,
ι
u
×
˙
v
=
ι
u
×
¨
ι
v
,
ι
u
/
˙
v
=
ι
u
/
¨
ι
v
,
u
≤
˙
v
⟺
ι
u
≤
¨
ι
v
.It turns out that ι
(
x
)
=
β
{
α
-
1
(
x
)
} for every x in A. Since, for example, u
+
˙
v
=
ι
-
1
{
ι
(
u
)
+
¨
ι
(
v
)
}, it should be clear that any statement in α-arithmetic can readily be transformed into a statement in β-arithmetic.


### 3.1.∗-Differentiation

If the following∗-limit in (8) exists, we denote it by f
∗
(
a
), call it the ∗-derivative of f at a, and say that f is ∗-differentiable at a: (8)
l
i
m
x
→
a
∗
f
x
-
¨
f
a
/
¨
ι
x
-
¨
ι
a
=
l
i
m
x
→
a
β
β
-
1
f
x
-
β
-
1
f
a
α
-
1
x
-
α
-
1
a
=
l
i
m
x
→
a
β
β
-
1
f
x
-
β
-
1
f
a
x
-
a
x
-
a
α
-
1
x
-
α
-
1
a
=
β
β
-
1
∘
f
′
a
α
-
1
′
a
.If the ∗-limit in (8) exists, then f
∗
(
a
) is necessarily in B. The ∗-derivative of f is the function that assigns to each number a in A the number f
∗
(
a
), if it exists. If, additionally, the second derivative of f at a exists, then by an easy substitution we obtain (9)
f
∗
∗
a
=
β
β
-
1
∘
f
′
a
α
-
1
′
a
′
1
α
-
1
′
a
=
β
β
-
1
∘
f
∗
a
′
α
-
1
′
a
.Repeating this procedure n times, we conclude that if its nth derivative at a exists, then f
∗
(
n
)
(
a
) exists and (10)
f
∗
n
a
=
β
β
-
1
∘
f
∗
n
-
1
a
′
α
-
1
′
a
,
n
=
1,2
,
…
.Let us formulate some useful results of differential calculus in terms of∗-derivative. They can be easily proved by application of the respective results of ∗-calculus to the function f.Theorem 1 (cf. [9] (mean value theorem)).
If the functionf is ∗-continuous on [
r
˙
,
s
˙
] and ∗-differentiable on (
r
˙
,
s
˙
), then there exists r
<
˙
c
<
˙
s such that (11)
f
∗
c
=
f
s
-
¨
f
r
ι
r
-
¨
ι
s
:
=
β
β
-
1
f
s
-
β
-
1
f
r
α
-
1
s
-
α
-
1
r
.


## 3.1.∗-Differentiation

If the following∗-limit in (8) exists, we denote it by f
∗
(
a
), call it the ∗-derivative of f at a, and say that f is ∗-differentiable at a: (8)
l
i
m
x
→
a
∗
f
x
-
¨
f
a
/
¨
ι
x
-
¨
ι
a
=
l
i
m
x
→
a
β
β
-
1
f
x
-
β
-
1
f
a
α
-
1
x
-
α
-
1
a
=
l
i
m
x
→
a
β
β
-
1
f
x
-
β
-
1
f
a
x
-
a
x
-
a
α
-
1
x
-
α
-
1
a
=
β
β
-
1
∘
f
′
a
α
-
1
′
a
.If the ∗-limit in (8) exists, then f
∗
(
a
) is necessarily in B. The ∗-derivative of f is the function that assigns to each number a in A the number f
∗
(
a
), if it exists. If, additionally, the second derivative of f at a exists, then by an easy substitution we obtain (9)
f
∗
∗
a
=
β
β
-
1
∘
f
′
a
α
-
1
′
a
′
1
α
-
1
′
a
=
β
β
-
1
∘
f
∗
a
′
α
-
1
′
a
.Repeating this procedure n times, we conclude that if its nth derivative at a exists, then f
∗
(
n
)
(
a
) exists and (10)
f
∗
n
a
=
β
β
-
1
∘
f
∗
n
-
1
a
′
α
-
1
′
a
,
n
=
1,2
,
…
.Let us formulate some useful results of differential calculus in terms of∗-derivative. They can be easily proved by application of the respective results of ∗-calculus to the function f.Theorem 1 (cf. [9] (mean value theorem)).
If the functionf is ∗-continuous on [
r
˙
,
s
˙
] and ∗-differentiable on (
r
˙
,
s
˙
), then there exists r
<
˙
c
<
˙
s such that (11)
f
∗
c
=
f
s
-
¨
f
r
ι
r
-
¨
ι
s
:
=
β
β
-
1
f
s
-
β
-
1
f
r
α
-
1
s
-
α
-
1
r
.


## 4. Non-Newtonian Differential Equations

We need some definitions and results from the theory of ordinary differential equations with respect to the∗-calculus before considering methods for approximating the solutions to initial value problems.Definition 2 (see [22]).
A functionf
(
t
,
y
) is said to satisfy the Lipschitz condition in the variable y on a set D
⊂
R
2 if a constant L
>
0 exists with (12)
f
t
,
y
1
-
f
t
,
y
2
≤
L
y
1
-
y
2whenever (
t
,
y
1
) and (
t
,
y
2
) are in D. The constant L is called a Lipschitz constant for f.Theorem 3.
Suppose thatD
=
{
(
t
,
y
)
∣
a
≤
˙
t
≤
˙
b
a
n
d
-
∞
<
¨
y
<
¨
∞
} and that f
(
t
,
y
) is ∗-continuous on D. If (
β
-
1
∘
f
) and β
′ satisfy a Lipschitz condition in the variable y on D and (
β
-
1
∘
f
) is bounded, then the initial value problem (13)
y
∗
t
=
f
t
,
y
t
,
a
≤
˙
t
≤
˙
b
,
y
a
=
λ
,has a unique solution y
(
t
) for a
≤
˙
t
≤
˙
b.Proof.
By using the inclusion (8) the ∗-differential equation in (13) can be transformed to the following ordinary differential equation: (14)
y
′
t
=
β
-
1
y
∗
t
β
′
y
t
α
′
t
,
α
′
t
≠
0
.Let us show that the function G
(
t
,
y
(
t
)
)
=
y
′
(
t
) in (14) satisfies the ordinary Lipschitz condition in y such that (15)
G
t
,
y
t
-
G
t
,
z
t
≤
β
′
y
t
α
′
t
β
-
1
f
t
,
y
t
-
β
′
y
t
α
′
t
β
-
1
f
t
,
z
t
+
β
′
y
t
α
′
t
β
-
1
f
t
,
z
t
-
β
′
z
t
α
′
t
β
-
1
f
t
,
z
t
≤
β
′
y
t
α
′
t
β
-
1
f
t
,
y
t
-
β
-
1
f
t
,
z
t
+
β
-
1
f
t
,
z
t
α
′
t
β
′
y
t
-
β
′
z
t
≤
β
′
y
t
α
′
t
L
1
y
t
-
z
t
+
β
-
1
f
t
,
z
t
α
′
t
L
2
y
t
-
z
t
≤
β
′
y
t
α
′
t
L
1
+
β
-
1
f
t
,
z
t
α
′
t
L
2
×
y
t
-
z
t
≤
L
y
t
-
z
t
,where L
1 and L
2 are Lipschitz constants of (
β
-
1
∘
f
) and β
′, respectively. We conclude that (13) has a unique solution. This step completes the proof.Theorem 4 (non-Newtonian Taylor’s theorem for one variable).
LetA be an open interval and let f
:
A
→
R be n
+
1 times ∗-differentiable on A. Then the function f can be represented in the neighborhood of a point a as a polynomial in t
-
˙
a, together with a remainder term, where c lies between t and a: (16)
f
t
=
β
︸
β
∑
k
=
0
n
‍
f
∗
k
a
×
¨
β
α
-
1
t
-
˙
a
k
k
!
+
¨
f
∗
n
+
1
c
×
¨
β
α
-
1
t
-
˙
a
n
+
1
n
+
1
!
.Bashirov et al. [9] stated non-Newtonian Taylor expansion of a function as (17)
f
t
=
∏
k
=
0
n
‍
f
∗
k
a
t
-
a
k
/
k
!
f
∗
n
+
1
c
t
-
a
n
+
1
/
n
+
1
!for α
=
I and β
=
e
x
p.Corollary 5 (∗-Taylor).
Iff
∗
(
n
)
(
t
) exists and is ∗-continuous for all n
∈
N at t
=
a and if l
i
m
∗
n
→
∞
R
n
(
t
)
=
β
(
0
), then, for all t
>
˙
a such that t
,
a
∈
R, (18)
f
t
=
l
i
m
n
→
∞
∗
β
︸
β
∑
k
=
0
n
‍
f
∗
k
a
×
¨
β
α
-
1
t
-
˙
a
k
k
!
.Following [22] we present non-Newtonian numerical methods which are a model of classical numerical methods, based on the concept of ∗-Taylor expansion.


### 4.1. Non-Newtonian Euler’s Method

Euler’s method is the most elementary approximation technique for solving initial value problems. Although it is seldom used in practice, the simplicity of its derivation can be used to illustrate the techniques involved in the construction of some of the more advanced techniques, without the cumbersome algebra that accompanies these constructions.The objective of Euler’s method is to obtain approximations to the initial value problem (13). A continuous approximation to the solution y
(
t
) will not be obtained; instead, approximations to y will be generated at various values, called mesh points, in the interval [
a
,
b
]. Once the approximate solution is obtained at the points, the approximate solution at other points in the interval can be found by interpolation.Selecting the mesh pointst
i
=
a
+
˙
i
×
˙
h for each i
=
0,1
,
2
,
…
,
N where the distance between the points h
=
t
i
+
1
-
˙
t
i is called the α-step size, we use Taylor’s theorem to derive Euler’s method. Suppose that y
(
t
), the unique solution to (13), has two continuous ∗-derivatives on [
a
˙
,
b
˙
], so that, for each i
=
0,1
,
2
,
…
,
N
-
1, (19)
y
t
i
+
1
=
y
t
i
+
¨
y
∗
t
i
×
¨
β
α
-
1
t
i
+
1
-
˙
t
i
+
¨
y
∗
∗
ξ
i
×
¨
β
α
-
1
t
i
+
1
-
˙
t
i
2
2
!
,for some number ξ
i
∈
t
˙
i
,
t
˙
i
+
1. Because h
=
t
i
+
1
-
˙
t
i, we have (20)
y
t
i
+
1
=
y
t
i
+
¨
f
t
i
,
y
t
i
×
¨
β
α
-
1
h
+
¨
y
∗
∗
ξ
i
×
¨
β
α
-
1
h
2
2
!
.Euler’s method constructs w
i
≈
y
(
t
i
), for each i
=
1,2
,
…
,
N
-
1, by deleting the remainder term. Thus, ∗-Euler method is (21)
w
0
=
λ
,
w
i
+
1
=
w
i
+
¨
β
α
-
1
h
×
¨
f
t
i
,
w
i
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
.If we choose the generators as α
=
e
x
p and β
=
e
x
p in (3), then (21) turns out to be (22)
w
0
=
λ
,
w
i
+
1
=
w
i
f
t
i
,
w
i
l
n
h
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
,and it is called bigeometric Euler method. Similarly we choose the generators as α
=
q and β
=
q in (4); then (21) turns out to be(23)
w
0
=
λ
,
w
i
+
1
=
w
i
p
+
h
f
t
i
,
w
i
p
1
/
p
,
p
∈
R
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
.and it is called q-Euler method.


### 4.2. Higher-Order∗-Taylor Methods

Since∗-Euler method was derived by using ∗-Taylor theorem with n
=
1 in (16) to approximate the solution of the differential equation, our first attempt to find methods for improving the convergence properties of difference methods is to extend this technique of derivation to larger values of n.Suppose the solutiony
(
t
) to the initial value problem (13) has n
+
1 continuous ∗-derivatives. If we expand the solution, y
(
t
), in terms of its nth ∗-Taylor polynomial about t
i and evaluate at t
i
+
1, we obtain (24)
y
t
i
+
1
=
y
t
i
+
¨
f
t
i
,
y
t
i
×
¨
β
α
-
1
h
+
¨
f
∗
t
i
,
y
t
i
×
¨
β
α
-
1
h
2
2
!
+
¨
⋯
+
¨
f
∗
n
ξ
i
,
y
ξ
i
×
¨
β
α
-
1
h
n
+
1
n
+
1
!for some ξ
i
∈
t
˙
i
,
t
˙
i
+
1. The difference-equation method corresponding to (24) is obtained by deleting the remainder term involving ξ
i. Taylor method constructs w
i
≈
y
(
t
i
), for each i
=
1,2
,
…
,
N
-
1, by deleting the remainder term. Thus ∗-Taylor method is (25)
w
0
=
λ
,
w
i
+
1
=
w
i
+
¨
T
n
t
i
,
w
i
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
,where (26)
T
n
t
i
,
w
i
=
β
α
-
1
h
×
¨
f
t
i
,
w
i
+
¨
β
α
-
1
h
2
2
!
×
¨
f
∗
t
i
,
w
i
+
¨
⋯
+
¨
β
α
-
1
h
n
n
!
×
¨
f
∗
n
-
1
t
i
,
w
i
.Therefore the bigeometric Taylor method can be obtained by choosing α
=
e
x
p and β
=
e
x
p in (25) as follows: (27)
w
0
=
λ
,
w
i
+
1
=
w
i
T
n
t
i
,
w
i
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
,where (28)
T
n
t
i
,
w
i
=
f
t
i
,
w
i
l
n
h
f
∗
t
i
,
w
i
l
n
h
2
/
2
⋯
f
∗
n
-
1
t
i
,
w
i
l
n
h
n
/
n
!
.Similarly q-Taylor method can be obtained by choosing α
=
q and β
=
q as follows: (29)
w
0
=
λ
,
w
i
+
1
=
w
i
p
+
T
n
t
i
,
w
i
p
1
/
p
,
p
∈
R
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
,where (30)
T
n
t
i
,
w
i
=
h
p
f
t
i
,
w
i
p
+
h
2
p
2
!
f
∗
t
i
,
w
i
p
v
v
v
v
+
⋯
+
h
n
p
n
!
f
∗
n
-
1
t
i
,
w
i
p
h
2
p
2
!
1
/
p
.


### 4.3. Non-Newtonian Runge-Kutta Method

Non-Newtonian Runge-Kutta (NNRK) methods have the high-order local truncation error of the∗-Taylor methods but eliminate the need to compute and evaluate the ∗-derivatives of f
(
t
,
y
). Before presenting the ideas behind their ∗-derivation, we need to consider Taylor theorem in two variables in the sense of non-Newtonian mean.Corollary 6 (mean value theorem).
Letf be a function of two variables x
=
x
(
t
) and y
=
y
(
t
) with continuous partial ∗-derivatives. Then there exist x
(
t
)
<
˙
c
1
<
˙
x
(
t
+
˙
Δ
t
) and y
(
t
)
<
¨
c
2
<
¨
y
(
t
+
¨
Δ
t
) such that (31)
f
x
∗
c
1
,
y
t
=
f
x
t
+
˙
Δ
t
,
y
t
-
¨
f
x
t
,
y
t
ι
x
t
+
˙
Δ
t
-
¨
ι
x
t
:
f
y
∗
x
t
,
c
2
=
f
x
t
,
y
t
+
¨
Δ
t
-
¨
f
x
t
,
y
t
ι
y
t
+
¨
Δ
t
-
¨
ι
y
t
:
.Theorem 7 (non-Newtonian chain rule for two variables).
Letf be a function of two variables x
=
x
(
t
) and y
=
y
(
t
) with continuous partial ∗-derivatives. If x and y are ∗-differentiable functions on (
a
˙
,
b
˙
) such that f
(
x
(
t
)
,
y
(
t
)
) is defined for every t
∈
(
a
˙
,
b
˙
), then (32)
d
∗
f
x
,
y
d
t
=
f
x
∗
x
,
y
×
¨
ι
∘
x
∗
t
+
¨
f
y
∗
x
,
y
×
¨
ι
∘
y
∗
t
.Proof.
LetΔ
x
=
x
(
t
+
˙
Δ
t
)
-
˙
x
(
t
) and Δ
y
=
y
(
t
+
¨
Δ
t
)
-
¨
y
(
t
). Then (33)
d
∗
f
x
,
y
d
t
=
l
i
m
Δ
t
→
0
∗
f
x
t
+
˙
Δ
t
,
y
t
+
˙
Δ
t
-
¨
f
x
,
y
ι
t
+
˙
Δ
t
-
¨
ι
t
:
=
l
i
m
Δ
t
→
0
∗
f
x
+
˙
Δ
x
,
y
+
¨
Δ
y
-
¨
f
x
,
y
ι
t
+
˙
Δ
t
-
¨
ι
t
:
=
l
i
m
∗
Δ
t
→
0
f
x
+
˙
Δ
x
,
y
-
¨
f
x
,
y
+
¨
f
x
+
˙
Δ
x
,
y
+
¨
Δ
y
-
¨
f
x
+
˙
Δ
x
,
y
ι
t
+
˙
Δ
t
-
¨
ι
t
:
=
l
i
m
Δ
t
→
0
∗
f
x
+
˙
Δ
x
,
y
-
¨
f
x
,
y
ι
t
+
˙
Δ
t
-
¨
ι
t
:
+
¨
l
i
m
Δ
t
→
0
∗
f
x
+
˙
Δ
x
,
y
+
¨
Δ
y
-
¨
f
x
+
˙
Δ
x
,
y
ι
t
+
˙
Δ
t
-
¨
ι
t
:.Combining this with relation (31) and (34)
l
i
m
Δ
t
→
0
∗
f
x
∗
c
1
,
y
=
f
x
∗
x
,
y
,
l
i
m
Δ
t
→
0
∗
f
y
∗
x
,
c
2
=
f
y
∗
x
,
y
,we get (35)
d
∗
f
x
,
y
d
t
=
l
i
m
Δ
t
→
0
∗
f
x
∗
c
1
,
y
×
¨
ι
x
t
+
˙
Δ
t
-
¨
ι
x
t
ι
t
+
˙
Δ
t
-
¨
ι
t
:
+
¨
l
i
m
Δ
t
→
0
∗
f
y
∗
x
,
c
2
×
¨
ι
y
t
+
˙
Δ
t
-
¨
ι
y
t
ι
t
+
˙
Δ
t
-
¨
ι
t
:
=
f
x
∗
x
,
y
×
¨
ι
∘
x
∗
t
+
¨
f
y
∗
x
,
y
×
¨
ι
∘
y
∗
t
.Theorem 8 (see [9] (non-Newtonian Taylor’s theorem for two variables)).
LetA be an open subset of R
2. Assume that the function f
:
A
→
R has all partial ∗-derivatives of order n
+
1 on A and let α and β be generator functions. Then, for every (
x
,
y
)
,
(
x
+
˙
h
,
y
+
¨
k
)
∈
A so that the line segment connecting these two points belongs to A, there exists a number θ
∈
(
0,1
) such that (36)
f
x
+
˙
h
,
y
+
¨
k
=
β
︸
β
∑
m
=
0
n
‍
β
︸
β
∑
i
=
0
m
‍
f
x
i
y
m
-
i
∗
m
x
,
y
×
¨
β
α
-
1
h
i
β
-
1
k
m
-
i
i
!
m
-
i
!
+
¨
β
︸
β
∑
i
=
0
n
+
1
‍
f
x
i
y
n
+
1
-
i
∗
n
+
1
x
+
˙
θ
h
,
y
+
¨
θ
k
×
¨
β
α
-
1
h
i
β
-
1
k
n
+
1
-
i
i
!
n
+
1
-
i
!
.Corollary 9.
Suppose thatf
(
t
,
y
) and all its partial ∗-derivatives of order less than or equal to n
+
1 are ∗-continuous on D
=
{
(
t
,
y
)
∣
a
≤
˙
t
≤
˙
b
,
c
≤
¨
t
≤
¨
d
}, and let (
t
0
,
y
0
)
∈
D. For every (
t
,
y
)
∈
D, there exists ξ between t and t
0 and μ between y and y
0 with (37)
f
t
,
y
=
P
n
t
,
y
+
¨
R
n
t
,
y
,where (38)
P
n
t
,
y
=
β
︸
β
∑
m
=
0
n
‍
β
︸
β
∑
i
=
0
m
‍
f
t
i
y
m
-
i
∗
m
t
0
,
y
0
×
¨
β
α
-
1
t
-
˙
t
0
i
β
-
1
y
-
¨
y
0
m
-
i
i
!
m
-
i
!
,
R
n
t
,
y
=
β
︸
β
∑
i
=
0
n
+
1
‍
f
t
i
y
n
+
1
-
i
∗
n
+
1
ξ
,
μ
×
¨
β
α
-
1
t
-
˙
t
0
i
β
-
1
y
-
¨
y
0
n
+
1
-
i
i
!
n
+
1
-
i
!
.The function P
n
(
t
,
y
) is called nth ∗-Taylor polynomial in two variables for the function f about (
t
0
,
y
0
), and R
n
(
t
,
y
) is the remainder term associated with P
n
(
t
,
y
).The first step in deriving a∗-Runge-Kutta method is to determine values for a
1
,
λ, and η with the property that a
1
×
¨
f
(
t
+
˙
λ
,
y
+
¨
η
) approximates (39)
T
2
t
,
y
=
β
α
-
1
h
×
¨
f
t
,
y
+
¨
β
α
-
1
h
2
2
×
¨
f
∗
t
,
y
.By taking into account the ∗-chain rule in Theorem 7, (40)
f
∗
t
,
y
=
f
t
∗
t
,
y
×
¨
ι
∘
t
∗
t
+
¨
f
y
∗
t
,
y
×
¨
ι
∘
y
∗
t
.We will from now on neglect the dependence of y on t when it appears as an argument to f. Therefore, the ∗-Taylor expansion becomes (41)
T
2
t
,
y
=
β
α
-
1
h
×
¨
f
t
,
y
+
¨
β
α
-
1
h
2
2
×
¨
f
t
∗
t
,
y
+
¨
f
y
∗
t
,
y
×
¨
ι
∘
y
∗
t
.Expanding f
(
t
+
˙
λ
,
y
+
¨
η
) in its ∗-Taylor polynomial of degree one about (
t
,
y
) in (36) gives (42)
a
1
×
¨
f
t
+
˙
λ
,
y
+
¨
η
=
a
1
×
¨
f
t
,
y
+
¨
a
1
×
¨
β
α
-
1
λ
×
¨
f
t
∗
t
,
y
+
¨
a
1
×
¨
β
β
-
1
η
×
¨
f
y
∗
t
,
y
+
¨
a
1
×
¨
R
1
t
+
˙
λ
,
y
+
¨
ηfor some ξ between t and t
+
˙
λ and μ between y and t
+
¨
η. Matching the coefficients of f and its ∗-derivatives gives the following three equations: (43)
a
1
=
β
α
-
1
h
,
a
1
×
¨
β
α
-
1
λ
=
β
α
-
1
h
2
2
,
a
1
×
¨
η
=
β
α
-
1
h
2
2
×
¨
ι
∘
y
∗
t
.The parameters a
1, λ, and η are therefore (44)
a
1
=
β
α
-
1
h
,
λ
=
α
α
-
1
h
2
,
η
=
β
α
-
1
h
β
-
1
ι
∘
y
∗
t
2
,so (45)
T
2
t
,
y
=
β
α
-
1
h
×
¨
f
t
+
˙
α
α
-
1
h
2
,
y
+
¨
β
α
-
1
h
β
-
1
ι
∘
y
∗
t
2
.Thus (46)
w
0
=
λ
,
w
i
+
1
=
w
i
+
¨
β
α
-
1
h
w
i
+
1
=
×
¨
f
t
i
+
˙
α
α
-
1
h
2
,
w
i
+
1
=
w
i
+
¨
β
α
-
1
h
β
-
1
ι
∘
y
∗
t
i
2for i
=
0,1
,
…
,
N
-
1. Only three parameters are present in a
1
×
¨
f
(
t
+
˙
λ
,
y
+
¨
η
) and all are needed in the match of T
(
2
). So a more complicated form is required to satisfy the conditions for any of the higher-order ∗-Taylor methods.Example 10 (see [8]).
As a simple and straight forward application of the second-order∗-Runge-Kutta method, we will consider the solution for the ordinary initial value problem (47)
y
′
t
=
1
-
1
t
,
y
1
=
1
.Obviously, the exact solution of (47) is y
(
t
)
=
t
-
l
n
t. We will show how second-order bigeometric Runge-Kutta (
B
G
R
K
2
) and q-Runge-Kutta (Q
R
K
2) methods work and will compare the results for appropriate step sizes.
Case 1. Before that, by taking into account the definition of the ∗-derivative given by (8), the corresponding bigeometric initial value problem for α
=
e
x
p, β
=
e
x
p is (48)
y
∗
t
=
e
x
p
t
y
′
t
y
t
=
e
x
p
t
-
1
y
,
y
1
=
1
.
Thus, we obtain from (46) (49)
y
t
h
=
y
t
f
t
h
,
y
∗
l
n
h
/
2
y
l
n
hwhich yields that (50)
w
0
=
λ
,
w
i
+
1
=
w
i
f
t
i
h
,
w
i
f
t
i
,
w
i
(
l
n
h
)
/
2
l
n
hfor i
=
0,1
,
…
,
N
-
1. Combining this with (48), we get the recurrence relation (51)
w
0
=
1
,
w
i
+
1
=
w
i
e
x
p
t
i
h
-
1
l
n
h
w
i
e
x
p
t
i
-
1
l
n
h
/
2
w
i
.
In Table1, the numerical results obtained by bigeometric non-Newtonian Runge-Kutta method and its comparison with the values of the exact solution with respect to the relative error are given. Due to the choice of generating function as exp, to be able to compute the values for x
>
x
0, one should choose step size h as h
>
1. As shown in the table in step 150 for h
=
2 we get the values near the point x
=
1
0
23 within very small relative error. Similarly, to compute the values for x
<
x
0, one should choose step size h as h
<
1. As it is seen from the table for h
=
0.99 in step 400 we get the values near the point x
=
0.01 and the relative error is quite small.
Furthermore, when we compare the results given in Table1 with the results of Rıza et al. (Tables 3-4) [8] where 4th-order bigeometric Runge-Kutta method and 4th-order ordinary Runge-Kutta method are compared with exact values, we see that the results given in Table 1 are fairly good. Therefore, we assert that for some family of problems the proposed non-Newtonian method gives better results.
Case 2. The initial value problem for α
=
q, β
=
q can be written as (52)
y
∗
t
=
y
1
-
1
/
p
t
-
1
1
/
p
t
,
y
1
=
1
.
Combining this with inclusion (46), (53)
y
t
p
+
h
p
1
/
p
=
t
p
+
h
p
2
1
/
p
,
y
p
+
h
p
f
p
t
,
y
2
1
/
p
p
y
p
+
h
p
f
v
v
v
v
×
t
p
+
h
p
2
1
/
p
,
y
p
+
h
p
f
p
t
,
y
2
1
/
p
p
1
/
p
.
Therefore, we derive that fori
=
0,1
,
…
,
N
-
1
(54)
w
0
=
1
,
w
i
+
1
=
w
i
p
+
1
2
h
p
w
i
p
+
h
p
w
i
p
-
1
t
i
-
1
2
t
i
p
1
-
1
/
p
∑
1
1
w
i
+
1
v
v
v
v
v
v
v
v
v
v
v
v
×
t
i
p
+
h
p
2
1
/
p
-
1
1
2
w
i
+
1
v
v
v
v
×
t
i
p
+
h
p
2
-
1
∑
1
1
1
/
p
.
In Table2, the numerical results obtained by Q
R
K
2 method and its comparison with the values of the exact solution with respect to the relative error are given. In this method since the step size depends on both h and p, for increasing values of p step size becomes smaller and smaller. Nevertheless, within close intervals, by taking different values of p and h, in particular (
p
=
2
,
h
=
0.01
) and (
p
=
10
,
h
=
0.5
), we obtained similar results. Moreover, when we compare the results in Table 1 for h
=
2 and the results for p
=
0.1 in Table 2, we see that B
G
R
K
2 method yields better results with respect to the results of Q
R
K
2.
Furthermore, when we compare the results given in Table2 with the results of Rıza et al. (Tables 3-4) [8] as before, we see that the results given in Table 2 are quite good.Table 1
Comparison of the results of the BGRK2 method with the exact values and its relative errors.
h
t
BGRK2
Exact
Relative error
2
t
1
=
2
1.067771
1.067639
1.2
×
1
0
-
4
t
10
=
2
10
28.54761
28.53426
4.7
×
1
0
-
4
t
100
=
2
100
0.1
×
1
0
16
0.1
×
1
0
16
0.1
×
1
0
-
13
t
150
=
2
150
0.3
×
1
0
23
0.3
×
1
0
23
0.2
×
1
0
-
18
0.99
t
1
=
0.99
1.000050
1.000050
0.4
×
1
0
-
7
t
100
=
0.36
1.371071
1.371065
0.4
×
1
0
-
5
t
250
=
0.08
2.593648
2.593642
0.2
×
1
0
-
5
t
400
=
0.01
4.038088
4.038084
0.9
×
1
0
-
6Table 2
Comparison of the results of the QRK2 method with the exact values and its relative errors.
p
h
t
QRK2
Exact
Relative error
2
0.5
t
5
=
1.5
1.095000
1.094534
0.000425
t
32
=
3.0
1.901872
1.901387
0.000255
t
140
=
6.0
4.208722
4.208240
0.000114
t
221
=
7.5
5.485578
5.485096
0.000087
2
0.01
t
1
=
1.000049
1.000000
1.000000
2.6
×
1
0
-
14
t
10
=
1.000499
1.000000
1.000000
2.5
×
1
0
-
13
t
21
=
1.001049
1.000000
1.000000
5.4
×
1
0
-
13
t
201
=
1.010000
1.000049
1.000049
0.5
×
1
0
-
11
10
0.5
t
5
=
1.000487
1.000000
1.000000
0.5
×
1
0
-
11
t
50
=
1.004778
1.000011
1.000011
0.5
×
1
0
-
10
t
150
=
1.013763
1.000093
1.000093
0.1
×
1
0
-
9
t
225
=
1.020061
1.000198
1.000198
0.1
×
1
0
-
9
0.1
0.01
t
12
=
0.2
×
1
0
10
0.3
×
1
0
10
0.2
×
1
0
10
0.494285
t
75
=
0.6
×
1
0
17
0.6
×
1
0
17
0.6
×
1
0
17
0.2
×
1
0
-
7
t
150
=
0.6
×
1
0
20
0.6
×
1
0
20
0.6
×
1
0
20
0.2
×
1
0
-
10
t
250
=
0.1
×
1
0
23
0.1
×
1
0
23
0.1
×
1
0
23
0.1
×
1
0
-
12Example 11 (see [23]).
Consider(55)
y
′
=
-
1
-
t
3
/
2
y
,
t
∈
0,1
,
y
0
=
1with solution y
(
t
)
=
e
x
p
{
2
/
5
[
(
1
-
t
)
5
/
2
-
1
]
}.
Case 1. The corresponding bigeometric initial value problem for α
=
e
x
p, β
=
e
x
p is (56)
y
∗
t
=
e
x
p
t
y
′
t
y
t
=
e
x
p
-
t
1
-
t
3
/
2
,
y
0
=
1
.
Combining this with (46), we get the recurrence relation (57)
w
0
=
1
,
w
i
+
1
=
w
i
e
x
p
-
t
i
h
1
-
t
i
h
3
/
2
l
n
h
.
In Table3, the numerical results obtained by bigeometric non-Newtonian Runge-Kutta method and its comparison with the values of the exact solution with respect to the relative error are given. The results were obtained for different values of h and we see that for smaller values of h the number of computations increases but relative error decreases. On the contrary, for bigger values of h, the number of computations decreases considerably but this leads to an increase in error.
Case 2. The initial value problem for α
=
q, β
=
q can be written as (58)
y
∗
x
=
-
y
p
-
1
1
-
t
3
/
2
y
t
p
-
1
1
/
p
,
y
1
=
1
.
Since the resulting recurrence relation is not defined at the pointt
=
0, we choose a starting value for t
>
0.
In Table4, the numerical results obtained by Q
R
K
2 method for p
=
2 and Q
R
K
2 method for p
=
1 which is the ordinary Runge-Kutta method are compared. As it is seen from the table the classical method yields better results. This shows that, in non-Newtonian approximation methods, the generating functions should be selected convenient to the problem under consideration. As shown in Tables 3 and 4 the B
G
R
K
2 method yields better results when it is compared with the results of the Q
R
K
2 method. This situation may be due to the exponential structure of the exact solution of the problem.Table 3
Comparison of the results of the BGRK2 method with the exact values and its relative errors.
h
t
BGRK2
Exact
Relative error
1
.
001
t
1
=
0.001
0.999000
0.999000
0.4
×
1
0
-
13
t
3000
=
0.020
0.980439
0.980439
0.7
×
1
0
-
9
t
6000
=
0.402
0.748647
0.748647
0.1
×
1
0
-
9
t
6911
=
0.999
0.670320
0.670320
0.1
×
1
0
-
8
1.1
t
50
=
0.001
0.998828
0.998827
0.4
×
1
0
-
6
t
75
=
0.012
0.987485
0.987481
0.4
×
1
0
-
5
t
115
=
0.575
0.702478
0.702521
0.6
×
1
0
-
4
t
120
=
0.927
0.670619
0.670105
0.1
×
1
0
-
3Table 4
Comparison of the results of the QRK2 method with the exact values and its relative errors.
p
h
t
QRK2
Exact
Relative error
2
0.05
t
1
=
0.050990
0.875026
0.874979
0.009667
t
100
=
0.500099
0.727159
0.719409
0.007750
t
250
=
0.790632
0.682999
0.675719
0.007279
t
399
=
0.998799
0.677541
0.670320
0.007221
1
0.05
t
1
=
0.150
0.907666
0.907666
0.4
×
1
0
-
4
t
8
=
0.500
0.719633
0.719434
0.1
×
1
0
-
3
t
14
=
0.800
0.675373
0.675133
0.2
×
1
0
-
3
t
18
=
1.000
0.670600
0.670320
0.2
×
1
0
-
3The areas of chemical engineering, nonlinear mechanics, biochemistry, and life sciences are sources of stiff problems. One possible application area of the present method can be the treatment of multiscale problems, that is, stiff ODEs, where different states reach equilibrium at different time scales.


## 4.1. Non-Newtonian Euler’s Method

Euler’s method is the most elementary approximation technique for solving initial value problems. Although it is seldom used in practice, the simplicity of its derivation can be used to illustrate the techniques involved in the construction of some of the more advanced techniques, without the cumbersome algebra that accompanies these constructions.The objective of Euler’s method is to obtain approximations to the initial value problem (13). A continuous approximation to the solution y
(
t
) will not be obtained; instead, approximations to y will be generated at various values, called mesh points, in the interval [
a
,
b
]. Once the approximate solution is obtained at the points, the approximate solution at other points in the interval can be found by interpolation.Selecting the mesh pointst
i
=
a
+
˙
i
×
˙
h for each i
=
0,1
,
2
,
…
,
N where the distance between the points h
=
t
i
+
1
-
˙
t
i is called the α-step size, we use Taylor’s theorem to derive Euler’s method. Suppose that y
(
t
), the unique solution to (13), has two continuous ∗-derivatives on [
a
˙
,
b
˙
], so that, for each i
=
0,1
,
2
,
…
,
N
-
1, (19)
y
t
i
+
1
=
y
t
i
+
¨
y
∗
t
i
×
¨
β
α
-
1
t
i
+
1
-
˙
t
i
+
¨
y
∗
∗
ξ
i
×
¨
β
α
-
1
t
i
+
1
-
˙
t
i
2
2
!
,for some number ξ
i
∈
t
˙
i
,
t
˙
i
+
1. Because h
=
t
i
+
1
-
˙
t
i, we have (20)
y
t
i
+
1
=
y
t
i
+
¨
f
t
i
,
y
t
i
×
¨
β
α
-
1
h
+
¨
y
∗
∗
ξ
i
×
¨
β
α
-
1
h
2
2
!
.Euler’s method constructs w
i
≈
y
(
t
i
), for each i
=
1,2
,
…
,
N
-
1, by deleting the remainder term. Thus, ∗-Euler method is (21)
w
0
=
λ
,
w
i
+
1
=
w
i
+
¨
β
α
-
1
h
×
¨
f
t
i
,
w
i
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
.If we choose the generators as α
=
e
x
p and β
=
e
x
p in (3), then (21) turns out to be (22)
w
0
=
λ
,
w
i
+
1
=
w
i
f
t
i
,
w
i
l
n
h
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
,and it is called bigeometric Euler method. Similarly we choose the generators as α
=
q and β
=
q in (4); then (21) turns out to be(23)
w
0
=
λ
,
w
i
+
1
=
w
i
p
+
h
f
t
i
,
w
i
p
1
/
p
,
p
∈
R
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
.and it is called q-Euler method.


## 4.2. Higher-Order∗-Taylor Methods

Since∗-Euler method was derived by using ∗-Taylor theorem with n
=
1 in (16) to approximate the solution of the differential equation, our first attempt to find methods for improving the convergence properties of difference methods is to extend this technique of derivation to larger values of n.Suppose the solutiony
(
t
) to the initial value problem (13) has n
+
1 continuous ∗-derivatives. If we expand the solution, y
(
t
), in terms of its nth ∗-Taylor polynomial about t
i and evaluate at t
i
+
1, we obtain (24)
y
t
i
+
1
=
y
t
i
+
¨
f
t
i
,
y
t
i
×
¨
β
α
-
1
h
+
¨
f
∗
t
i
,
y
t
i
×
¨
β
α
-
1
h
2
2
!
+
¨
⋯
+
¨
f
∗
n
ξ
i
,
y
ξ
i
×
¨
β
α
-
1
h
n
+
1
n
+
1
!for some ξ
i
∈
t
˙
i
,
t
˙
i
+
1. The difference-equation method corresponding to (24) is obtained by deleting the remainder term involving ξ
i. Taylor method constructs w
i
≈
y
(
t
i
), for each i
=
1,2
,
…
,
N
-
1, by deleting the remainder term. Thus ∗-Taylor method is (25)
w
0
=
λ
,
w
i
+
1
=
w
i
+
¨
T
n
t
i
,
w
i
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
,where (26)
T
n
t
i
,
w
i
=
β
α
-
1
h
×
¨
f
t
i
,
w
i
+
¨
β
α
-
1
h
2
2
!
×
¨
f
∗
t
i
,
w
i
+
¨
⋯
+
¨
β
α
-
1
h
n
n
!
×
¨
f
∗
n
-
1
t
i
,
w
i
.Therefore the bigeometric Taylor method can be obtained by choosing α
=
e
x
p and β
=
e
x
p in (25) as follows: (27)
w
0
=
λ
,
w
i
+
1
=
w
i
T
n
t
i
,
w
i
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
,where (28)
T
n
t
i
,
w
i
=
f
t
i
,
w
i
l
n
h
f
∗
t
i
,
w
i
l
n
h
2
/
2
⋯
f
∗
n
-
1
t
i
,
w
i
l
n
h
n
/
n
!
.Similarly q-Taylor method can be obtained by choosing α
=
q and β
=
q as follows: (29)
w
0
=
λ
,
w
i
+
1
=
w
i
p
+
T
n
t
i
,
w
i
p
1
/
p
,
p
∈
R
f
o
r
e
a
c
h
i
=
0,1
,
2
,
…
,
N
-
1
,where (30)
T
n
t
i
,
w
i
=
h
p
f
t
i
,
w
i
p
+
h
2
p
2
!
f
∗
t
i
,
w
i
p
v
v
v
v
+
⋯
+
h
n
p
n
!
f
∗
n
-
1
t
i
,
w
i
p
h
2
p
2
!
1
/
p
.


## 4.3. Non-Newtonian Runge-Kutta Method

Non-Newtonian Runge-Kutta (NNRK) methods have the high-order local truncation error of the∗-Taylor methods but eliminate the need to compute and evaluate the ∗-derivatives of f
(
t
,
y
). Before presenting the ideas behind their ∗-derivation, we need to consider Taylor theorem in two variables in the sense of non-Newtonian mean.Corollary 6 (mean value theorem).
Letf be a function of two variables x
=
x
(
t
) and y
=
y
(
t
) with continuous partial ∗-derivatives. Then there exist x
(
t
)
<
˙
c
1
<
˙
x
(
t
+
˙
Δ
t
) and y
(
t
)
<
¨
c
2
<
¨
y
(
t
+
¨
Δ
t
) such that (31)
f
x
∗
c
1
,
y
t
=
f
x
t
+
˙
Δ
t
,
y
t
-
¨
f
x
t
,
y
t
ι
x
t
+
˙
Δ
t
-
¨
ι
x
t
:
f
y
∗
x
t
,
c
2
=
f
x
t
,
y
t
+
¨
Δ
t
-
¨
f
x
t
,
y
t
ι
y
t
+
¨
Δ
t
-
¨
ι
y
t
:
.Theorem 7 (non-Newtonian chain rule for two variables).
Letf be a function of two variables x
=
x
(
t
) and y
=
y
(
t
) with continuous partial ∗-derivatives. If x and y are ∗-differentiable functions on (
a
˙
,
b
˙
) such that f
(
x
(
t
)
,
y
(
t
)
) is defined for every t
∈
(
a
˙
,
b
˙
), then (32)
d
∗
f
x
,
y
d
t
=
f
x
∗
x
,
y
×
¨
ι
∘
x
∗
t
+
¨
f
y
∗
x
,
y
×
¨
ι
∘
y
∗
t
.Proof.
LetΔ
x
=
x
(
t
+
˙
Δ
t
)
-
˙
x
(
t
) and Δ
y
=
y
(
t
+
¨
Δ
t
)
-
¨
y
(
t
). Then (33)
d
∗
f
x
,
y
d
t
=
l
i
m
Δ
t
→
0
∗
f
x
t
+
˙
Δ
t
,
y
t
+
˙
Δ
t
-
¨
f
x
,
y
ι
t
+
˙
Δ
t
-
¨
ι
t
:
=
l
i
m
Δ
t
→
0
∗
f
x
+
˙
Δ
x
,
y
+
¨
Δ
y
-
¨
f
x
,
y
ι
t
+
˙
Δ
t
-
¨
ι
t
:
=
l
i
m
∗
Δ
t
→
0
f
x
+
˙
Δ
x
,
y
-
¨
f
x
,
y
+
¨
f
x
+
˙
Δ
x
,
y
+
¨
Δ
y
-
¨
f
x
+
˙
Δ
x
,
y
ι
t
+
˙
Δ
t
-
¨
ι
t
:
=
l
i
m
Δ
t
→
0
∗
f
x
+
˙
Δ
x
,
y
-
¨
f
x
,
y
ι
t
+
˙
Δ
t
-
¨
ι
t
:
+
¨
l
i
m
Δ
t
→
0
∗
f
x
+
˙
Δ
x
,
y
+
¨
Δ
y
-
¨
f
x
+
˙
Δ
x
,
y
ι
t
+
˙
Δ
t
-
¨
ι
t
:.Combining this with relation (31) and (34)
l
i
m
Δ
t
→
0
∗
f
x
∗
c
1
,
y
=
f
x
∗
x
,
y
,
l
i
m
Δ
t
→
0
∗
f
y
∗
x
,
c
2
=
f
y
∗
x
,
y
,we get (35)
d
∗
f
x
,
y
d
t
=
l
i
m
Δ
t
→
0
∗
f
x
∗
c
1
,
y
×
¨
ι
x
t
+
˙
Δ
t
-
¨
ι
x
t
ι
t
+
˙
Δ
t
-
¨
ι
t
:
+
¨
l
i
m
Δ
t
→
0
∗
f
y
∗
x
,
c
2
×
¨
ι
y
t
+
˙
Δ
t
-
¨
ι
y
t
ι
t
+
˙
Δ
t
-
¨
ι
t
:
=
f
x
∗
x
,
y
×
¨
ι
∘
x
∗
t
+
¨
f
y
∗
x
,
y
×
¨
ι
∘
y
∗
t
.Theorem 8 (see [9] (non-Newtonian Taylor’s theorem for two variables)).
LetA be an open subset of R
2. Assume that the function f
:
A
→
R has all partial ∗-derivatives of order n
+
1 on A and let α and β be generator functions. Then, for every (
x
,
y
)
,
(
x
+
˙
h
,
y
+
¨
k
)
∈
A so that the line segment connecting these two points belongs to A, there exists a number θ
∈
(
0,1
) such that (36)
f
x
+
˙
h
,
y
+
¨
k
=
β
︸
β
∑
m
=
0
n
‍
β
︸
β
∑
i
=
0
m
‍
f
x
i
y
m
-
i
∗
m
x
,
y
×
¨
β
α
-
1
h
i
β
-
1
k
m
-
i
i
!
m
-
i
!
+
¨
β
︸
β
∑
i
=
0
n
+
1
‍
f
x
i
y
n
+
1
-
i
∗
n
+
1
x
+
˙
θ
h
,
y
+
¨
θ
k
×
¨
β
α
-
1
h
i
β
-
1
k
n
+
1
-
i
i
!
n
+
1
-
i
!
.Corollary 9.
Suppose thatf
(
t
,
y
) and all its partial ∗-derivatives of order less than or equal to n
+
1 are ∗-continuous on D
=
{
(
t
,
y
)
∣
a
≤
˙
t
≤
˙
b
,
c
≤
¨
t
≤
¨
d
}, and let (
t
0
,
y
0
)
∈
D. For every (
t
,
y
)
∈
D, there exists ξ between t and t
0 and μ between y and y
0 with (37)
f
t
,
y
=
P
n
t
,
y
+
¨
R
n
t
,
y
,where (38)
P
n
t
,
y
=
β
︸
β
∑
m
=
0
n
‍
β
︸
β
∑
i
=
0
m
‍
f
t
i
y
m
-
i
∗
m
t
0
,
y
0
×
¨
β
α
-
1
t
-
˙
t
0
i
β
-
1
y
-
¨
y
0
m
-
i
i
!
m
-
i
!
,
R
n
t
,
y
=
β
︸
β
∑
i
=
0
n
+
1
‍
f
t
i
y
n
+
1
-
i
∗
n
+
1
ξ
,
μ
×
¨
β
α
-
1
t
-
˙
t
0
i
β
-
1
y
-
¨
y
0
n
+
1
-
i
i
!
n
+
1
-
i
!
.The function P
n
(
t
,
y
) is called nth ∗-Taylor polynomial in two variables for the function f about (
t
0
,
y
0
), and R
n
(
t
,
y
) is the remainder term associated with P
n
(
t
,
y
).The first step in deriving a∗-Runge-Kutta method is to determine values for a
1
,
λ, and η with the property that a
1
×
¨
f
(
t
+
˙
λ
,
y
+
¨
η
) approximates (39)
T
2
t
,
y
=
β
α
-
1
h
×
¨
f
t
,
y
+
¨
β
α
-
1
h
2
2
×
¨
f
∗
t
,
y
.By taking into account the ∗-chain rule in Theorem 7, (40)
f
∗
t
,
y
=
f
t
∗
t
,
y
×
¨
ι
∘
t
∗
t
+
¨
f
y
∗
t
,
y
×
¨
ι
∘
y
∗
t
.We will from now on neglect the dependence of y on t when it appears as an argument to f. Therefore, the ∗-Taylor expansion becomes (41)
T
2
t
,
y
=
β
α
-
1
h
×
¨
f
t
,
y
+
¨
β
α
-
1
h
2
2
×
¨
f
t
∗
t
,
y
+
¨
f
y
∗
t
,
y
×
¨
ι
∘
y
∗
t
.Expanding f
(
t
+
˙
λ
,
y
+
¨
η
) in its ∗-Taylor polynomial of degree one about (
t
,
y
) in (36) gives (42)
a
1
×
¨
f
t
+
˙
λ
,
y
+
¨
η
=
a
1
×
¨
f
t
,
y
+
¨
a
1
×
¨
β
α
-
1
λ
×
¨
f
t
∗
t
,
y
+
¨
a
1
×
¨
β
β
-
1
η
×
¨
f
y
∗
t
,
y
+
¨
a
1
×
¨
R
1
t
+
˙
λ
,
y
+
¨
ηfor some ξ between t and t
+
˙
λ and μ between y and t
+
¨
η. Matching the coefficients of f and its ∗-derivatives gives the following three equations: (43)
a
1
=
β
α
-
1
h
,
a
1
×
¨
β
α
-
1
λ
=
β
α
-
1
h
2
2
,
a
1
×
¨
η
=
β
α
-
1
h
2
2
×
¨
ι
∘
y
∗
t
.The parameters a
1, λ, and η are therefore (44)
a
1
=
β
α
-
1
h
,
λ
=
α
α
-
1
h
2
,
η
=
β
α
-
1
h
β
-
1
ι
∘
y
∗
t
2
,so (45)
T
2
t
,
y
=
β
α
-
1
h
×
¨
f
t
+
˙
α
α
-
1
h
2
,
y
+
¨
β
α
-
1
h
β
-
1
ι
∘
y
∗
t
2
.Thus (46)
w
0
=
λ
,
w
i
+
1
=
w
i
+
¨
β
α
-
1
h
w
i
+
1
=
×
¨
f
t
i
+
˙
α
α
-
1
h
2
,
w
i
+
1
=
w
i
+
¨
β
α
-
1
h
β
-
1
ι
∘
y
∗
t
i
2for i
=
0,1
,
…
,
N
-
1. Only three parameters are present in a
1
×
¨
f
(
t
+
˙
λ
,
y
+
¨
η
) and all are needed in the match of T
(
2
). So a more complicated form is required to satisfy the conditions for any of the higher-order ∗-Taylor methods.Example 10 (see [8]).
As a simple and straight forward application of the second-order∗-Runge-Kutta method, we will consider the solution for the ordinary initial value problem (47)
y
′
t
=
1
-
1
t
,
y
1
=
1
.Obviously, the exact solution of (47) is y
(
t
)
=
t
-
l
n
t. We will show how second-order bigeometric Runge-Kutta (
B
G
R
K
2
) and q-Runge-Kutta (Q
R
K
2) methods work and will compare the results for appropriate step sizes.
Case 1. Before that, by taking into account the definition of the ∗-derivative given by (8), the corresponding bigeometric initial value problem for α
=
e
x
p, β
=
e
x
p is (48)
y
∗
t
=
e
x
p
t
y
′
t
y
t
=
e
x
p
t
-
1
y
,
y
1
=
1
.
Thus, we obtain from (46) (49)
y
t
h
=
y
t
f
t
h
,
y
∗
l
n
h
/
2
y
l
n
hwhich yields that (50)
w
0
=
λ
,
w
i
+
1
=
w
i
f
t
i
h
,
w
i
f
t
i
,
w
i
(
l
n
h
)
/
2
l
n
hfor i
=
0,1
,
…
,
N
-
1. Combining this with (48), we get the recurrence relation (51)
w
0
=
1
,
w
i
+
1
=
w
i
e
x
p
t
i
h
-
1
l
n
h
w
i
e
x
p
t
i
-
1
l
n
h
/
2
w
i
.
In Table1, the numerical results obtained by bigeometric non-Newtonian Runge-Kutta method and its comparison with the values of the exact solution with respect to the relative error are given. Due to the choice of generating function as exp, to be able to compute the values for x
>
x
0, one should choose step size h as h
>
1. As shown in the table in step 150 for h
=
2 we get the values near the point x
=
1
0
23 within very small relative error. Similarly, to compute the values for x
<
x
0, one should choose step size h as h
<
1. As it is seen from the table for h
=
0.99 in step 400 we get the values near the point x
=
0.01 and the relative error is quite small.
Furthermore, when we compare the results given in Table1 with the results of Rıza et al. (Tables 3-4) [8] where 4th-order bigeometric Runge-Kutta method and 4th-order ordinary Runge-Kutta method are compared with exact values, we see that the results given in Table 1 are fairly good. Therefore, we assert that for some family of problems the proposed non-Newtonian method gives better results.
Case 2. The initial value problem for α
=
q, β
=
q can be written as (52)
y
∗
t
=
y
1
-
1
/
p
t
-
1
1
/
p
t
,
y
1
=
1
.
Combining this with inclusion (46), (53)
y
t
p
+
h
p
1
/
p
=
t
p
+
h
p
2
1
/
p
,
y
p
+
h
p
f
p
t
,
y
2
1
/
p
p
y
p
+
h
p
f
v
v
v
v
×
t
p
+
h
p
2
1
/
p
,
y
p
+
h
p
f
p
t
,
y
2
1
/
p
p
1
/
p
.
Therefore, we derive that fori
=
0,1
,
…
,
N
-
1
(54)
w
0
=
1
,
w
i
+
1
=
w
i
p
+
1
2
h
p
w
i
p
+
h
p
w
i
p
-
1
t
i
-
1
2
t
i
p
1
-
1
/
p
∑
1
1
w
i
+
1
v
v
v
v
v
v
v
v
v
v
v
v
×
t
i
p
+
h
p
2
1
/
p
-
1
1
2
w
i
+
1
v
v
v
v
×
t
i
p
+
h
p
2
-
1
∑
1
1
1
/
p
.
In Table2, the numerical results obtained by Q
R
K
2 method and its comparison with the values of the exact solution with respect to the relative error are given. In this method since the step size depends on both h and p, for increasing values of p step size becomes smaller and smaller. Nevertheless, within close intervals, by taking different values of p and h, in particular (
p
=
2
,
h
=
0.01
) and (
p
=
10
,
h
=
0.5
), we obtained similar results. Moreover, when we compare the results in Table 1 for h
=
2 and the results for p
=
0.1 in Table 2, we see that B
G
R
K
2 method yields better results with respect to the results of Q
R
K
2.
Furthermore, when we compare the results given in Table2 with the results of Rıza et al. (Tables 3-4) [8] as before, we see that the results given in Table 2 are quite good.Table 1
Comparison of the results of the BGRK2 method with the exact values and its relative errors.
h
t
BGRK2
Exact
Relative error
2
t
1
=
2
1.067771
1.067639
1.2
×
1
0
-
4
t
10
=
2
10
28.54761
28.53426
4.7
×
1
0
-
4
t
100
=
2
100
0.1
×
1
0
16
0.1
×
1
0
16
0.1
×
1
0
-
13
t
150
=
2
150
0.3
×
1
0
23
0.3
×
1
0
23
0.2
×
1
0
-
18
0.99
t
1
=
0.99
1.000050
1.000050
0.4
×
1
0
-
7
t
100
=
0.36
1.371071
1.371065
0.4
×
1
0
-
5
t
250
=
0.08
2.593648
2.593642
0.2
×
1
0
-
5
t
400
=
0.01
4.038088
4.038084
0.9
×
1
0
-
6Table 2
Comparison of the results of the QRK2 method with the exact values and its relative errors.
p
h
t
QRK2
Exact
Relative error
2
0.5
t
5
=
1.5
1.095000
1.094534
0.000425
t
32
=
3.0
1.901872
1.901387
0.000255
t
140
=
6.0
4.208722
4.208240
0.000114
t
221
=
7.5
5.485578
5.485096
0.000087
2
0.01
t
1
=
1.000049
1.000000
1.000000
2.6
×
1
0
-
14
t
10
=
1.000499
1.000000
1.000000
2.5
×
1
0
-
13
t
21
=
1.001049
1.000000
1.000000
5.4
×
1
0
-
13
t
201
=
1.010000
1.000049
1.000049
0.5
×
1
0
-
11
10
0.5
t
5
=
1.000487
1.000000
1.000000
0.5
×
1
0
-
11
t
50
=
1.004778
1.000011
1.000011
0.5
×
1
0
-
10
t
150
=
1.013763
1.000093
1.000093
0.1
×
1
0
-
9
t
225
=
1.020061
1.000198
1.000198
0.1
×
1
0
-
9
0.1
0.01
t
12
=
0.2
×
1
0
10
0.3
×
1
0
10
0.2
×
1
0
10
0.494285
t
75
=
0.6
×
1
0
17
0.6
×
1
0
17
0.6
×
1
0
17
0.2
×
1
0
-
7
t
150
=
0.6
×
1
0
20
0.6
×
1
0
20
0.6
×
1
0
20
0.2
×
1
0
-
10
t
250
=
0.1
×
1
0
23
0.1
×
1
0
23
0.1
×
1
0
23
0.1
×
1
0
-
12Example 11 (see [23]).
Consider(55)
y
′
=
-
1
-
t
3
/
2
y
,
t
∈
0,1
,
y
0
=
1with solution y
(
t
)
=
e
x
p
{
2
/
5
[
(
1
-
t
)
5
/
2
-
1
]
}.
Case 1. The corresponding bigeometric initial value problem for α
=
e
x
p, β
=
e
x
p is (56)
y
∗
t
=
e
x
p
t
y
′
t
y
t
=
e
x
p
-
t
1
-
t
3
/
2
,
y
0
=
1
.
Combining this with (46), we get the recurrence relation (57)
w
0
=
1
,
w
i
+
1
=
w
i
e
x
p
-
t
i
h
1
-
t
i
h
3
/
2
l
n
h
.
In Table3, the numerical results obtained by bigeometric non-Newtonian Runge-Kutta method and its comparison with the values of the exact solution with respect to the relative error are given. The results were obtained for different values of h and we see that for smaller values of h the number of computations increases but relative error decreases. On the contrary, for bigger values of h, the number of computations decreases considerably but this leads to an increase in error.
Case 2. The initial value problem for α
=
q, β
=
q can be written as (58)
y
∗
x
=
-
y
p
-
1
1
-
t
3
/
2
y
t
p
-
1
1
/
p
,
y
1
=
1
.
Since the resulting recurrence relation is not defined at the pointt
=
0, we choose a starting value for t
>
0.
In Table4, the numerical results obtained by Q
R
K
2 method for p
=
2 and Q
R
K
2 method for p
=
1 which is the ordinary Runge-Kutta method are compared. As it is seen from the table the classical method yields better results. This shows that, in non-Newtonian approximation methods, the generating functions should be selected convenient to the problem under consideration. As shown in Tables 3 and 4 the B
G
R
K
2 method yields better results when it is compared with the results of the Q
R
K
2 method. This situation may be due to the exponential structure of the exact solution of the problem.Table 3
Comparison of the results of the BGRK2 method with the exact values and its relative errors.
h
t
BGRK2
Exact
Relative error
1
.
001
t
1
=
0.001
0.999000
0.999000
0.4
×
1
0
-
13
t
3000
=
0.020
0.980439
0.980439
0.7
×
1
0
-
9
t
6000
=
0.402
0.748647
0.748647
0.1
×
1
0
-
9
t
6911
=
0.999
0.670320
0.670320
0.1
×
1
0
-
8
1.1
t
50
=
0.001
0.998828
0.998827
0.4
×
1
0
-
6
t
75
=
0.012
0.987485
0.987481
0.4
×
1
0
-
5
t
115
=
0.575
0.702478
0.702521
0.6
×
1
0
-
4
t
120
=
0.927
0.670619
0.670105
0.1
×
1
0
-
3Table 4
Comparison of the results of the QRK2 method with the exact values and its relative errors.
p
h
t
QRK2
Exact
Relative error
2
0.05
t
1
=
0.050990
0.875026
0.874979
0.009667
t
100
=
0.500099
0.727159
0.719409
0.007750
t
250
=
0.790632
0.682999
0.675719
0.007279
t
399
=
0.998799
0.677541
0.670320
0.007221
1
0.05
t
1
=
0.150
0.907666
0.907666
0.4
×
1
0
-
4
t
8
=
0.500
0.719633
0.719434
0.1
×
1
0
-
3
t
14
=
0.800
0.675373
0.675133
0.2
×
1
0
-
3
t
18
=
1.000
0.670600
0.670320
0.2
×
1
0
-
3The areas of chemical engineering, nonlinear mechanics, biochemistry, and life sciences are sources of stiff problems. One possible application area of the present method can be the treatment of multiscale problems, that is, stiff ODEs, where different states reach equilibrium at different time scales.


## 5. Concluding Remarks

Although there are many excellent ways to present the principle of the classical calculus, the novel presentations probably lead most naturally to the development of the non-Newtonian calculi in general. Every property in classical calculus has an analogue in non-Newtonian calculus which is a methodology that allows one to have a different look at problems which can be investigated via calculus.In this paper it is shown that, depending on the structure of the problem, even theN
N
R
K
2 method gives favorable results with respect to classical Runge-Kutta method. Based on the selection of different generating functions, the approximations within the same method vary according to the problem under consideration. The numerical results are obtained for narrower and larger intervals in the same table and it is acquired as an advantage of the non-Newtonian method. In particular, for large intervals the selection of “exp” function and for small intervals the selection of “q” function as a generator is seen to be more appropriate. Therefore, depending on the selection of the appropriate generating functions, closer approximations can be made in macro- and microsystems for real-life situations.We see that theN
N
R
K
2 method is an applicable tool for the solution of initial value problems and that, because of the nature of the ∗-derivative, the results may be superior to the ones obtained from the ordinary derivative. The family of problems where the higher-order NNRK method gives better results to the N
N
R
K
2 method has to be studied explicitly as a future project for the numerical solutions of partial differential equations in [24–26].



---

*Source: 594685-2015-02-05.xml*