Chapter 4. Differential equation and applications. (c) 2001 Akio Wakabayashi

Chapter 4.. .Differential equations and applications.

1. Differential equations of the first order.
1.1 The method of solution by this operational calculus.
. . The differential equations of the first order can generally be expressed in the following equation, where f(y) depends on the function exciting the differential equation but the expression is omitted.

________(1.1) Supposing the values of the derivative y' (t) and the function f(y) are y'₀, y'₁, y'₂, ……… and f₀, f₁, f₂, ……… at the equidistant points of t_k=t₀+kh, there exist the relations that,
Eq1_1a

These relations can be expressed in the vector of this operational calculus as follows.
Eq1_2

______(1.2) This is the differential equation which denotes the differential equation in Eq. (1.1). The vector F is of the vector-valued function expressed in F(Y ). Using the differential operator, the differential equation may be expressed in,
Eq1_3

______(1.3). . The solution of this differential equation can be expressed in the result of multiplying the both sides by the integral operator hD⁻¹ or h(ΔD)⁻¹, that is,
Eq1_4

______(1.4)
Eq1_5

______(1.5) If the vector-valued function F is independent of the solution Y, the solution is obtained only with the product of the operator and F. However, the accuracy of the solution depends on the value of the equidistant interval. In order to solve the differential equation with ill condition, the iterative method must be used by transforming it into the differential equation whose right side function depends on the solution y(t). The examples are mentioned in Chapter 2 Section 2.7.
. . When the right side function f depends on the solution y(t), the integral operation in Eq. (1.4) has the following relation for any Y and Z.
Eq1_6

In the same way, the integral operation in Eq. (1.5) also has the following result. Eq1_7

______(1.7) Hence, if the condition,
Eq1_8

________(1.8) is satisfied, the mappings denoted by the product of the integral operator and the operand are both contraction and the equations in Eq. (1.4) and (1.5) have a unique solution by fixed point theorem. The solution can be obtained by the iterative method. Denoting the starting value and the iterative results by Y₀ and Y₁, Y₂, Y₃, ………, the iteration is carried out as follows.
Eq1_8a

Because of the condition (1.7) and (1.8), these iteration converges on the solution Y, that is,
Eq1_9

________(1.9) The starting value Y₀ is the vector whose elements are all zero except the first element is y₀. The vector F(Y₀) also is the vector whose elements are all zero except the first element is f(y₀). These elements being all zero may be omitted because it has no meaning to carry out the calculation of these elements. Accordingly, the dimensions of the vector Y and F grow by iteration as represented below.
. . If the initial value y₀ is zero and the value f(y₀) also is zero and the function f(y) has not the variable t except of y(t), the iterative method gives the solution which is identical with zero. If there exists a solution which is not identical with zero, the first iterative procedure must slightly be modified to obtain the solution. If one of the above conditions is not satisfied, the solution being not identical with zero is obtained.
Eq1_9a

Δy₀ is corrected and Δ²y₀ is predicted. Eq1_9c

4)__In the same way, the correction with prediction is iterated until Y_n−1 is obtained, where n is predetermined.
Eq1_9d

Δy₀, Δ²y₀, ………, Δⁿ⁻¹y₀ are corrected and Δⁿy₀ is predicted. Therefore,
Eq1_9e

6)__All the elements of Y_n are corrected by the integral operator which has the rectangular matrix of order (n×(n+1)).
Eq1_10

___(1.10)
. . If the norm of Y_n1−Y_n is less than the value ε, which is larger than the smallest value of the last significant decimal place of the numerical solution, the result Y_n1 is the convergent solution. If the condition is not satisfied, there is the need to iterate the correction expressed in Eq. (1.10) once more by use of Y_n1 instead of Y_n. Denoting the result by Y_n2, if the norm of Y_n2−Y_n1 is less than the value ε, the result Y_n2 is the convergent solution. If the condition is not satisfied the result is not convergent yet. In the case, there is a need to increase the value n or to decrease the equidistant interval h if the condition (1.8) is not satisfied sufficiently or if the check is not made. In order to save the time spent on calculation, it is better to make the equidistant interval half.
. . In case of n=2 and n=4, the author showed the programs which does not use the vector calculus but transforms above vector equations into the expressions using the values y_k, f_k of the functions at equidistant points. However, it is better to use the vector calculus when the n value is larger than three, because the subroutine computing the product of the first row of the integral matrix and the operand vector can be used in order to compute all the products of the other rows and the operand vector. It makes the program simple and makes the parallel processing easy. The program and application are represented in Chapter 2 Section 5.

1.2 The solution by symbolical formula manipulation.
. . Denoting the equidistant interval by the symbol h, the elements of vectors are expressed in use of the symbol h. Using the integral matrix whose elements are the symbolical expressions in Eq. (1.10), the symbolical formula manipulation can obtain the solution expressed in the symbolical formula, which is Taylor's expansion with remainder, where the remainder has some errors.
. . The differential equation y' =−ty, y(0)=1 has the solution expressed in Taylor's expansion with remainder as,

______(1.11) because the differential coefficients of successive order at t=0 are as follows.
Eq1_11a

. . When we solve this differential equation by the method in the preceding section, the symbolical formula manipulation becomes complex if f(y) is transformed into the vector F by calculating f(y)=−ty(t) after transforming the vector Y into the function y(t). If f(y) is the function of y only, the vector F can be obtained by the vectorial calculus from the vector Y immediately. However, in the case that f(y) includes explicitly the independent variable t as this differential equation, the independent variable t also must be transformed into the vector T because it is also the function of t. The vector T is expressed in the vector T of Eq. (1.12) because t are the equidistant points from t=0 to some point with interval h. By this transform, the variable t of function y(t), in which t is not expressed explicitly also is transformed into the vector T. Accordingly, both vector Y and F becomes the function of vector T. However, the vector F is not changed by this transform, so whichever the vector F is the function of t or the function of T, it has no influence on the differential operator and the integral operator. When the vector F is the function of independent variable t, the author named it the vector-valued function. When the vector F is the function of vector T, the author names it the vectorial function. By this transform, the numerical calculus also can calculate the function F(Y) on the right side immediately from the vector Y without transforming the vector Y into the values of function y(t).
. . Using the vectorial function, this differential equation is expressed in the vector equations as follows.
Eq1_12

_____(1.12) Premultiplying the both sides of the differential equation by the integral operator hD⁻¹ and expanding the vector T into the matrix T, the solution is expressed in,
Eq1_13

____(1.13). . Starting the iterative method by use of the starting value y₀1, the solution is obtained as follows, if the n value is three.
1)__ Eq1_13a

______The error is the remainder y' (ξ)h.
Eq1_13b

4)__Correcting Y₃ by use of the integral matrix of order (3×4),
Eq1_13d

5)__Correcting Y₄ by use of the integral matrix of order (3×4) once more,
Eq1_13e

. . The solution expressed in the symbolical formula is the identical part between the results obtained by the twice correction. More than twice corrections using the integral matrix of order (3×4) do not yield the correct coefficient of the term with h⁶ and yield many terms which have h⁸, h⁹, ……… and whose coefficients are very large. These terms result in irregular numerical solution if the interval 3h does not sufficiently satisfy the condition by which Taylor's expansion converges.However, it must be taken care that all the terms that are different from the higher order terms of Taylor's expansion are not always errors in practically numerical solutions. In order that the numerical solution is exact value, the remainder term is required, adding to the terms that are the same as Taylor's expansion with h⁴ and lower order. These higher order terms may be a part of the remainder term, so these are not always numerical error. The examples are shown in the numerical examples of the extrapolation method in the section 1.4.
. . In the case of Taylor's expansion in Eq. (1.11), calculating the successive differences of y(t) at t=0, h, 2h, 3h, the following results are obtained and comparing each element of the vector solved above with the following results, the terms on and after h⁴ are not equal. However, the term h⁴ of the step 4) is equal to the term h⁴ of the step 5), so the term h⁴ is convergent. Furthermore, calculating the solutions y₁, y₂, y₃ by use of the terms up to the term h⁴, these are equal to Taylor's expansion Eq. (1.11) up to the term h⁴. This shows that the term h⁴ need not be equal to Taylor's expansion and that it need be convergent. Accordingly, un-transforming the vector into the function, the deciding of convergence can be carried out by use of the elements of vector. In numerical solution too, it is able to use the same way.
Eq1_13f1

. . In numerical calculation, the results of the step 3), 4) and 5) give the data by which it is able to decide whether the numerical solution converges or not. If there is no significant difference between the results of the step 3) and 4), the difference between the coefficients of the terms with h⁴ is not significant and the solution has been convergent. If there is some significant difference between them, the difference between the coefficients of the terms with h⁴ is significant. In the case, if there is no significant difference between the results of the step 4) and 5), the difference between the coefficients of the terms with h⁶ is not significant and the solution has been convergent. If there is some significant difference between them, the difference between the coefficients of the terms with h⁶ is significant and the solution is not convergent yet. In the case, if the interval 4h sufficiently satisfies the condition required by Taylor's expansion, the n value may be increased by one. If the interval 4h does not sufficiently satisfy the condition, the equidistant interval h must be decreased. These are the reason why the author's operational calculus uses the procedure mentioned in the preceding section.
. . If the equidistant interval has been decreased the iterative method must be restarted from the step 1). If the n value has been increased it must be restarted from the step 1), or it must be continued from the step 4) by appending the fourth row of the integral matrix D⁻¹. The result of the step 5) must not be used, because it has the terms which have h⁶ and h⁸ and whose coefficients are irregular. Accordingly, in case of n=4, the steps on and after 4) are modified as,
4)__Appending the fourth row of the integral matrix, Δ⁴y₀ is appended to ΔY₄ as,
Eq1_13g

5)__Correcting Y₄ by use of the integral matrix of order (4×5),
Eq1_13h

6)__Correcting Y₅ by use of the integral matrix of order (4×5) once more,
Eq1_13j

. . The predicting step 4) by use of the integral matrix of order four gives the result that only y₂ has the correct terms up to h⁶. However, the correcting step 5) by use of the integral matrix of order (4×5) gives the result that only y₄ has the correct terms up to h⁶. The step 6) correcting once more by use of the integral matrix of order (4×5) gives the result that all the terms up to h⁶ are convergent but the coefficients of the terms of h⁶ are not correct except y₄ has the correct terms of h⁶ and h⁸. In order to obtain all the correct terms of h⁶, we must carried out the predicting step 5) by use of the integral matrix of order five and twice correcting steps by use of the integral matrix of order (5×6).
. . These are the reason why the programs in Chapter 2 judge the convergence of the numerical solutions by the convergence of y_n instead of the norm in case of using the integral matrix of order (n×(n+1)) and why they start to solve the next n intervals by use of y_n as the initial value. The progress by one interval as usual predictor-corrector methods loses the accuracy and stability of the solution.
. . Milne's predictor is of the same as the calculation of y₄ in the above step 4) but its term of h⁶ is not correct. His corrector cannot make the term correct because it is Simpson's 1/3 rule, whose precision is worse than Milne's predictor. His predictor is obtained by use of the integral matrix of order four in step 5) of the preceding section as follows.

_____(1.14) Substituting the values given by the differential equation y' =−ty into Eq. (1.14) and giving y₀=1 and three accurate values y₁, y₂, y₃ expressed in Eq. (1.11),
Eq1_15

___(1.15) Comparing this result with that of author's operational calculus shown above step 4), it has the remainders of y₁, y₂ and y₃ together with the remainder of y₄ as the cause of error. This shows that if y₁, y₂ and y₃ have some errors yielded by calculation, these errors cause the error of y₄, too.
. . Substituting the values given by the differential equation into Milne's corrector and giving two accurate value y₂, y₃ expressed in Eq. (1.11) and the predicted value y₄, where the error except of the term of h⁶ is denoted by ER,
Eq1_16

____(1.16) The coefficient of h⁶ is not correct but it has been convergent because it is determined by the coefficient of h⁶ of y₂ and the coefficients of h⁴ of y₂, y₃, y₄ and these are correct and convergent value already. Hence, Milne's corrector has no effect of the iterations except it yields many cause of error when the differential equation has ill condition.
. . Usually Milne's method is regarded as the method integrating the interval h from t₃=t₀+3h to t₄=t₀+4h by use of the results y₀, y₁, y₂, y₃. However, it is the integral of the four intervals from t₀ to t₄ as shown in Eq. (1.15) and Eq. (1.16). In order to make the series convergent, it is required that the four intervals 4h sufficiently satisfy the condition which must make Taylor's expansion convergent. On the other hand, Runge-Kutta's method integrates the interval h by use of the half interval as shown below and requires that only the interval h must sufficiently satisfy the condition.
Eq1_16a

Accordingly, usual comparison of Milne's method and Runge-Kutta's method is not correct. The comparison must be made by the condition that Milne's method integrates the interval h by use of the four equidistant interval of h/4 in case of comparison with usual Runge-Kutta's method and that Runge-Kutta's method integrates the interval 4h by use of the half interval 2h in case of comparison with usual Milne's method.
. . Table 1.1 shows the comparison of results in numerical calculation, where the initial value is y(0)=10. The column named Milne 1 shows the results of Milne's method using the corrector only once. The column named Milne 2 shows the results in case of using the corrector N times more as usual method, where the N values are shown in brackets and denote the convergence of solution except N=21 denotes no convergence. The one is more accurate and stable than the other after t=5. The iterative use of the corrector has almost no effect except of bud consequences after t=5.
. . The column named R.-K. 1 shows the results of Runge-Kutta's method integrating four intervals, which are the integral range of Milne's predictor. The results are slightly worse than Milne 1 and have no meaning after t=4. The column named R.-K. 2 shows the results of Runge-Kutta's method integrating two intervals, which are the integral range of Milne's corrector. The results are slightly better than Milne 1 around t=0 but they become worse than Milne 1 after t=1.6. The column named R.-K. 3 shows the results of Runge-Kutta's method integrating one interval, which is the integral range used by usual comparison of the two methods. The results are better than Milne 1 but the comparison has no meaning because the integral range is a half of Milne's corrector.

Table 1.1 Comparison of Milne's method with Runge-Kutta's method.
t Milne 1 [N] Milne 2 True 10ⁿ R.-K. 1 R.-K. 2 R.-K.3

-0.2
-0.1
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2.0
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
3.0
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
4.0
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
5.0
5.1
5.2
5.3
5.4
5.5
5.6
5.7
5.8
5.9
6.0 *
*
*
*
9.801985
9.559971
9.231156
8.824958
8.352687
7.827025
7.261467
6.669740
6.065277
5.460712
4.867491
4.295543
3.753084
3.246501
2.780355
2.357447
1.978980
1.644742
1.353355
1.102511
8.892256
7.100641
5.613601
4.393810
3.404867
2.612241
1.984204
1.492146
1.110959
8.189017
5.976289
4.317867
3.088738
2.187330
1.533693
1.064518
7.316307
4.976955
3.352948
2.235188
1.476127
9.641744
6.242481
3.994103
2.535664
1.588600
9.891442
6.065259
3.704269
2.222832
1.330509
7.824053
4.573649
2.649826
1.498494
8.658934
4.659715
2.722967
1.388222 *
*
*
*
[1] 9.801985
[2] 9.559971
[2] 9.231157
[2] 8.824961
[2] 8.352690
[2] 7.827032
[2] 7.261475
[2] 6.669752
[2] 6.065289
[2] 5.460727
[21] 4.867506
[2] 4.295559
[2] 3.753097
[3] 3.246514
[2] 2.780364
[3] 2.357455
[3] 1.978982
[3] 1.644743
[3] 1.353352
[3] 1.102507
[4] 8.892171
[4] 7.100566
[3] 5.613492
[4] 4.393728
[3] 3.404758
[3] 2.612172
[4] 1.984113
[3] 1.492105
[3] 1.110893
[4] 8.188920
[21] 5.975855
[5] 4.318048
[5] 3.088463
[21] 2.187724
[5] 1.533482
[5] 1.065061
[5] 7.313908
[21] 4.983469
[4] 3.349478
[7] 2.242724
[6] 1.470929
[8] 9.730303
[7] 6.166814
[8] 4.102402
[8] 2.428532
[9] 1.726982
[8] 8.395819
[9] 7.900475
[9] 1.620723
[10] 4.726138
[10]-1.587930
[11] 4.269(-5)
[21]-3.669(-5)
[11] 5.203(-5)
[11]-5.751(-5)
[11] 7.179(-5)
[11]-8.492(-5)
[12] 1.035(-4)
[21]-1.249(-4) 9.80198673
9.95012479
1.00000000
9.95012479
9.80198673
9.55997478
9.23116344
8.82496903
8.35270199
7.82704545
7.26149030
6.66976789
6.06530660
5.46074412
4.86752228
4.29557318
3.75311111
3.24652467
2.78037290
2.35746057
1.97898674
1.64474464
1.35335283
1.10250492
8.89216081
7.10053615
5.61347500
4.39369336
3.40474421
2.61214065
1.98410974
1.49207819
1.11089965
8.18869738
5.97602198
4.31784069
3.08871441
2.18749112
1.53380989
1.06476605
7.31802551
4.97955236
3.35462628
2.23745881
1.47748183
9.65933345
6.25214775
4.00652974
2.54193577
1.59667624
9.92949522
6.11356511
3.72665317
2.24905706
1.34381028
7.94938558
4.65571332
2.69957850
1.54975396
8.80816483
4.95639984
2.76124090
1.52299797 +0
+0
+1
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
-1
-1
-1
-1
-1
-1
-1
-1
-1
-2
-2
-2
-2
-2
-2
-2
-3
-3
-3
-3
-3
-4
-4
-4
-4
-4
-5
-5
-5
-5
-5
-6
-6
-6
-6
-7
-7
-7
-7 *
*
*
*
*
*
9.231147
*
*
*
7.261394
*
*
*
4.867729
*
*
*
2.782125
*
*
*
1.357819
*
*
*
5.681144
*
*
*
2.055314
*
*
*
6.536229
*
*
*
1.879862
*
*
*
5.107530
*
*
*
1.390(-3)
*
*
*
4.047(-4)
*
*
*
1.345(-4)
*
*
*
5.365(-5)
*
*
*
2.659(-5) *
*
*
*
9.801987
*
9.231163
*
8.352701
*
7.261490
*
6.065314
*
4.867548
*
3.753172
*
2.780485
*
1.979161
*
1.353590
*
8.895055
*
5.616689
*
3.408020
*
1.987196
*
1.113606
*
5.998208
*
3.105769
*
1.546143
*
7.402126
*
3.408833
*
1.510565
*
6.443620
*
2.647193
*
1.047986
*
4.000694
*
1.473920
*
5.245412
*
1.805212
*
6.015464
*
1.943720 *
*
*
9.950125
9.801987
9.559975
9.231163
8.824968
8.352701
7.827044
7.261489
6.669767
6.065306
5.460744
4.867524
4.295576
3.753114
3.246529
2.780379
2.357469
1.978997
1.644756
1.353366
1.102520
8.892320
7.100702
5.613648
4.393867
3.404918
2.612308
1.984269
1.492228
1.111037
8.189947
5.977133
4.318812
3.089554
2.188204
1.534407
1.065258
7.322023
4.982759
3.357159
2.239432
1.479001
9.670855
6.260771
4.012900
2.546580
1.600024
9.953292
6.130269
3.738234
2.256986
1.349177
7.985244
4.679385
2.715018
1.559705
8.871569
4.996307
2.786066
1.538262

. . In case of the integral of two intervals, the author's method gives following results expressed in symbolical formula by modifying the step 3) and 4) into correcting procedures using the integral matrix of order (2×3). The expression y₂ of the step 4) is of the same as Runge-Kutta's result integrating two intervals 2h by use of the interval h as a half interval. Accordingly, the comparison of numerical results must be made with R.-K. 2.
Eq1_16b
Eq1_16c
. . The program equivalent to this method is represented in Chapter 2 Section 1.3, where the vector calculus is not used. Rewriting the subroutines on and after the line 56 for solving above differential equation from zero to six by use of 30 equidistant intervals, it gives the result y=1.522434×10⁻⁷ at t=6. The result is better than R.-K. 3, which integrates a half of integral interval of the author's method. This shows that Runge-Kutta's method yields more errors in numerical calculation.
. . The author's method also is superior to usual methods by the reason that the results y₂ of the step 2), 3) and 4) are the data by which we can determine the convergence of the solution y₂ or the suitableness of integral interval and that we can obtain more accurate result by varying the integral interval. The author's variable interval method gives the result y=2.005017×10⁻³⁶ at t=13.0 as shown in Table 2.7 of Chapter 2 Section 2.4. The true value of the solution is 2.00500878……×10⁻³⁶.

1.3 Lipschitz condition.
. . The value of integral operator is constant value, so the left side of Eq. (1.7) becomes,

Therefore, Eq. (1.7) is_______________ Eq1_17

____________(1.17)
Eq. (1.8) requires that there exists a constant k such that |F' (Y)|≤k is satisfied. Hence the vector-valued function F(Y) or the function f(y) must satisfy Lipschitz condition. However, Lipschitz constant is max|F' (ξ)| or max|f ' (ξ)| strictly by the procedure obtaining Eq. (1.7), so the function f(y) may have a point where the value f ' (y) is infinity if there is a mean value max|f ' (ξ)| between y and z.
. . It is said that the differential equation in Eq. (1.18) cannot be solved by numerical calculus because the function f(y), that is, y' (y) does not satisfy Lipschitz condition as Eqs. (1.19).
Eq1_18

__________(1.18)
Eq1_19

____(1.19) However, the function f(y) satisfies Lipschitz condition because even lf one of y and z is zero there exists the mean value f' (ξ) as follows.
Eq1_19a

. . The numerical solution of this differential equation becomes y(t)≡0 because Runge-Kutta's results are all zero by the following results and because Milne's method and other predictor-corrector methods use those.
Eq1_19b

. . In quadrature, we obtain the solution by supposing y(t)>0 and dividing the both sides of the differential equation by the right side. The numerical methods must be carried out on the supposition too. The trapezoidal method is able to solve the differential equation by use of iterative method on the supposition that an arbitrary value ε is the approximate solution at t=h.
Eq1_19c

This result is the fixed point of the equation as to y₁ for y₁>0 as follows. Eq1_20

______(1.20) Accordingly, the equation expressed in the integral for Eq. (1.18), that is,
Eq1_21

________(1.21) has the fixed point when the integral is approximated by the trapezoidal method, even if the h value is larger than one. Hence, the function f(y) satisfies Lipschitz condition. The initial value at t=0 has no influence on Lipschitz condition, because the equation in Eq. (1.21) has the initial value y₀ in the both cases of approximate y and z. The initial values also have no influence on Eq. (1.17) because of y₀=z₀.
. . The equation in Eq. (1.21) satisfies the condition Eq. (1.8) required in order to have the fixed point, independently of the h value. Supposing y>ξ>z>0,
Eq1_21a

Accordingly, it is independent of Lipschitz condition that usual methods cannot solve the differential equation in Eq. (1.18) and it is only dependent on f₀=0 at y₀=0.
. . The author's method is able to solve the differential equation by use of the fixed point in Eq. (1.20) as the integrand of the step 2), skipping the step 1).
Eq1_21b

Usually, it is said that the solution of the differential equation in Eq. (1.18) is,
Eq1_22

for arbitrary c.______(1.22) In the case, there exists the hidden condition that the solution is not always zero when the t value is equal to or larger than an arbitrary constant c. The numerical method also is able to obtain the solution in Eq. (1.22), if it uses the hidden condition.
. . The above example is the case that the fixed point of the trapezoidal method is the accurate solution. Even if the fixed point is the approximate solution, the author's method is able to obtain the numerically correct result. In case of the differential equation in Eq. (1.12), the fixed point of the trapezoidal method is expressed in the series whose terms on and after h⁴ are not correct. Eq1_22a

Predicting ΔY by use of the integral matrix of order (2×2) and correcting it twice by use of the integral matrix of order (2×3),
Eq1_22b

The result y₂ is of the same as the result of the step 4) in the last of the previous section, if the terms on and after h⁸ are negligible numerically.
. . The numerical methods cannot solve the differential equation in Eq. (1.23), because the value f(0) is infinity. The reason is independent of Lipschitz condition because there exists the mean value f ' (ξ) as shown in Eq. (1.24), if yz≠0.
Eq1_23

________(1.23)
Eq1_24

_____(1.24) In quadrature, we obtain the solution by multiplying the both sides of the differential equation by y, supposing y(t)≠0 and then we suppose it satisfies the initial condition y(0)=0. The author's method also is able to use the technique, that is, it can solve the differential equation by giving the initial value at t=ε, which is smaller than the h value.
. . The trapezoidal method also is able to solve the differential equation if the initial value y₀ is given at t=h when h is the equidistant interval. Supposing the value y₁ is an arbitrary value at t=2h and iterating the trapezoidal method in Eq. (1.25), the solutions are given as shown in column A of Table 1.2, when the equidistant interval is h=0.1. Eq1_25

______(1.25)

Table 1.2__y₀=y(0.1), h=0.1
t	A	B	True value
0.1 0.2 0.3 0.4 0.5 : 1.0 2.0 3.0 4.0 5.0 : : 10.0 : : 100.0	0.447213 0.637454 0.779994 0.899673 1.00500 : 1.41819 2.00297 2.45196 2.83058 3.16421 : : 4.47352 : : 14.1426	0.450000 0.639319 0.781506 0.900981 1.00617 : 1.41902 2.00355 2.45243 2.83100 3.16459 : : 4.47378 : : 14.1427	(Initial value) 0.6324555 0.7745966 0.8944271 1.0000000 : 1.4142135 2.0000000 2.4494897 2.8284271 3.1622776 : : 4.4721359 : : 14.142135

. .
These are equivalent to the fixed points of the following equations as to y_k.
Eq1_26

_____(1.26). . The trapezoidal method loses many significant figures of the solution y₁ at t=2h, whereas the accurate initial value is given at t=h. However, it gets the lost significant figures according as the t value increases. It means that the trapezoidal method does not require the initial value with many significant figures but two significant figures as shown in column B of Table 1.2. Accordingly, if we can estimate the approximate value with the two significant figures at t=h by use of the initial value y=0 at t=0, we can solve the differential equation by the initial value.
. . The mean value of f₀ and f₁ of Eq. (1.25) is very different from the mean value of integral around t=0.1. The exact mean value f is given by the reciprocal of f₀ and f₁ as follows. Changing the independent variable into y by the reciprocal of the both sides of the differential equation in Eq. (1.23),
Eq1_26a

__________(1.27) Using the mean value f instead of the mean value of Eq. (1.25), the exact solution is obtained by iterative method. Supposing the solution y₁ at t=2h is equal to the initial value y₀ at t=h, the iterative solution converges to the exact value as follows.

0.4472136, 0.6708204, 0.6260990, 0.6335526, 0.6322675, 0.6324878, 0.6324500
0.6324565, 0.6324554, 0.6324556, 0.6324555

. . The mean value f in Eq. (1.27) is held even if the initial value y₀ is zero at t=0. Accordingly, supposing the solution y₁ is an arbitrary value c at t=h, the mean value f of integral is 2/c and y₀=y₁−f₁h=c−2h/c=0 must be satisfied. Therefore, the solution y₁ is the fixed point of the equation as to c in Eq. (1.28). However, the iterative method obtains,
Eq1_27a

The sequence of y₁ is bistable and not convergent. This iteration uses y₀ as fixed value, so the cause that the sequence does not converge is irrelevant to Lipschitz condition. It is due to the reason that this iteration uses the mean value of integral in Eq. (1.27) supposing y₀=0 whereas y₀=y₁−f₁h does not become zero. The iteration using the mean value given by this non-zero y₀ and Eq. (1.27) does not converge by the reason that the point whose value is fixed does not exist and that Lipschitz condition is not satisfied at t=0. Above sequence y₁ is the result of right iteration and the equation in Eq. (1.28) shows that though the bistable values are different formally these are equal value essentially. Accordingly, even if the bistable values are numeric, the exact solution of y₁ is obtained by supposing the product of the bistable values becomes the square of y₁. The difference of the bistable values varies according to the arbitrary value c given by supposition. If the two values of the bistable are equal the difference becomes zero. In the case, the c value is the exact solution y₁ at t=h and it is obtained by the equation in Eq. (1.28). In numerical calculus, supposing the mean value of the bistable values as the corrected value of y₁, the iteration converges to the exact solution. By use of these ideas, the differential equation in Eq. (1.23) can be solved numerically from t=0.
Eq1_28

________(1.28) If the function f(y) is the reciprocal of yⁿ,
Eq1_29

. . In general, these differential equations may be solved as the solutions t_k at the equidistant points y_k=kh by the author's methods. However, if we would like to obtain the solutions y_k at the equidistant points t_k=kh, we may solve it by the method using variable interval as follows.

1) We obtain the solution t=δ at y=ε, which is very small.

2) Changing the independent variable into t, we obtain the solution y₁ by the integral from t=δ to t=h of f(y).

3) We obtain the solution y_k at t=kh, (k>1) with the initial value y₁ at t=h.

. . In case of the differential equation in Eq. (1.23), supposing ε is 0.1, the half interval is y=0.05 and the step 1) gives,

The step 2) gives the solution y₁ at t=h=0.1 as in Table 1.3 by use of Program 2.4 in Chapter 2 Section 2.7, where the subroutines on and after Line 56 are,

56 F=1/Y(1):RETURN
58 XB=0.005:XE=0.1:N=1:ND=1:Y(1)=0.1:EE=16:RETURN
60 Y#=SQR(X#):PRINT USING"True=##.########^^^^";Y#:RETURN

Line 58 instructs to solve one interval from 0.005 to 0.1 by subdividing it by 16. The interval is subdivided by 64 automatically because the solution does not converge at the first equidistant point with the interval 0.095/16 and 0.095/32. When the solution has arrived at the tenth interval, the intervals after the tenth interval are changed into the twice interval, that is, 0.095/32. When the solution has arrived at t=0.1, the interval is 0.095/8. The true value of the solution at t=0.1 is different from the exact value at the eighth decimal place but it is due to the error of converting the value 0.1 into the binary system with single precision.
. . The step 3) gives the solutions y_k+1 at t=0.1+kh for k=IB as in Table 1.4, where Line 58 is,
58 XB=0.1:XE=10:N=99:ND=1:Y(1)=0.4472137:EE=8:RETURN The line 58 instructs to solve the 99 equidistant intervals from 0.1 to 10 with the length 0.1 by subdividing one interval by 8. The subdivided interval becomes twice interval after t=0.175, 0.35, 0.7 and the interval becomes 0.1 after t=0.7. The true value of the solution at t=0.5 is not one but it is due to the rounded error of the t value calculated. The t value is 0.499999970 and it is equivalent to 0.5 in BASIC with single precision.

Table 1.3__The solution from t=0.005 to 0.1.

EB EE t (10ⁿ) y (10ⁿ) True (10ⁿ)

1
2
3
4
5
6
7
8
9
10
6
7
8
9
10
6
7
8
9
10
11
12
7
8 64
64
64
64
64
64
64
64
64
64
32
32
32
32
32
16
16
16
16
16
16
16
8
8 6.48438(-3)
7.96875(-3)
9.45313(-3)
0.0109375
0.0124219
0.0139062
0.0153906
0.016875
0.0183594
0.0198437
0.0228125
0.0257813
0.02875
0.0317187
0.0346875
0.040625
0.0465625
0.0525
0.0584375
0.064375
0.0703125
0.07625
0.088125
0.1 1.138805(-1)
1.262439(-1)
1.375001(-1)
1.479020(-1)
1.576190(-1)
1.667709(-1)
1.754459(-1)
1.837118(-1)
1.916214(-1)
1.992173(-1)
2.136001(-1)
2.270738(-1)
2.397916(-1)
2.518681(-1)
2.633914(-1)
2.850439(-1)
3.051640(-1)
3.240371(-1)
3.418699(-1)
3.588176(-1)
3.750001(-1)
3.905126(-1)
4.198215(-1)
4.472137(-1) 1.13880419(-1)
1.26243811(-1)
1.37500002(-1)
1.47901993(-1)
1.57619004(-1)
1.66770799(-1)
1.75445862(-1)
1.83711734(-1)
1.91621377(-1)
1.99217217(-1)
2.13600096(-1)
2.27073783(-1)
2.39791578(-1)
2.51868019(-1)
2.63391346(-1)
2.85043851(-1)
3.05163891(-1)
3.24037030(-1)
3.41869859(-1)
3.58817497(-1)
3.75000000(-1)
3.90512488(-1)
4.19821403(-1)
4.47213599(-1)

Table 1.4__The solution from t=0.1 to 10.

IB EB EE t y (10ⁿ) True (10ⁿ)

1 1
2
3
4
5
6
4 8
8
8
8
8
8
4 0.1125
0.125
0.1375
0.15
0.1625
0.175
0.2 4.743421(-1)
5.000004(-1)
5.244048(-1)
5.477229(-1)
5.700880(-1)
5.916082(-1)
6.324558(-1) 4.74341658(-1)
5.00000000(-1)
5.24404430(-1)
5.47722568(-1)
5.70087702(-1)
5.91607973(-1)
6.32455537(-1)

2 1
2
3
4 4
4
4
4 0.225
0.25
0.275
0.3 6.708207(-1)
7.071070(-1)
7.416201(-1)
7.745969(-1) 6.70820407(-1)
7.07106781(-1)
7.41619857(-1)
7.74596646(-1)

3 1
2
2 4
4
2 0.325
0.35
0.4 8.062260(-1)
8.366602(-1)
8.944274(-1) 8.06225760(-1)
8.36660019(-1)
8.94427164(-1)

4 1
2 2
2 0.45
0.5 9.486835(-1)
1.000000 9.48683285(-1)
9.99999970(-1)

5 1
2 2
2 0.55
0.6 1.048809
1.095445 1.04880880
1.09544508

6 1
2 2
2 0.65
0.7 1.140176
1.183216 1.14017540
1.18321595

7
8
9
19
29
39
49
59
69
79
89
99 1
1
1
1
1
1
1
1
1
1
1
1 1
1
1
1
1
1
1
1
1
1
1
1 0.8
0.9
1
2
3
4
5
6
7
8
9
10 1.264911
1.341641
1.414214
2.000000
2.449490
2.828427
3.162278
3.464102
3.741658
4.000001
4.242641
4.472137 1.26491107
1.34164077
1.41421356
1.99999994
2.44948965
2.82842704
3.16227751
3.46410148
3.74165726
3.99999988
4.24264069
4.47213595

1.4 The solution by interpolation and extrapolation.
. . The author's methods using the integral matrix of order n require n+1 equidistant points. Usually, these points are obtained by dividing the each integral interval by n as mentioned in previous section, when the integral interval is given by dividing the given integral range by the number of points at which the solutions should be obtained. Hence the author names them the solution by interpolation. It is similar to Runge-Kutta's method, which uses the predicted result at the half point of the integral interval except the way of using the result is different from the author's method using 3 points.
. . The solution by interpolation is able to solve the differential equation by use of the variable interval into which every interval subdivided and by use of the integral matrix of variable order. However, it cannot integrate tabulated functions because it cannot obtain the values at equidistant points dividing the tabulated interval by n. In this case, the author's methods must use n tabulated interval to obtain the required n+1 equidistant points. The author names these the solution by extrapolation because these use the points t=t₀+kh for k>1 in order to correct the solution at t=t₀+h.
. . The solution by extrapolation is able to use the integral matrix of a fixed order or variable order. The both methods are able to use two ways advancing the lower limit of integral. One advances to t₀+nh if the integral matrix is of the order (n×(n+1)), when the differential equation has well condition such that the solution converges at t₀+nh. Another advances to t₀+kh for 1≤k≤n, when the solutions are convergent at the successive k points. This way is very effective when the integral matrix is of variable order.
. . Table 1.5 shows the solution of the differential equation y' =−ty, y(0)=10 by the extrapolation using the equidistant interval h=0.1 and the integral matrix of variable order n≤11. The column IB is the equidistant points and the column E is the order of the integral matrix, that is, the corrector is of order ((E−1)×E). When IB is one, the solution at t=0.1 converges by use of the four equidistant points t=0, 0.1, 0.2, 0.3 and the solutions at t=0.2, 0.3 are not convergent, so the IB value is increased by 1. When IB is two, all the solutions converge at t=0.2, 0.3, 0.4, 0.5 by use of the solution at t=0.1 as the initial value, so the IB value is increased by four. When IB is equal to or larger than 28, the solution converges only at the point denoted by the IB value. When IB is equal to or larger than 52, the order of the integral matrix does not vary by the restriction given. The restriction loses the accuracy of solution. However, the loss is very slight, compared with Milne's method and Runge-Kutta's method shown in Table 1.1.

Table 1.5 . .The solution by extrapolation of y' =−ty, y(0)=10.
IB	E	t	y (10ⁿ)	True (10ⁿ)	. .	IB	E	t	y (10ⁿ)	True (10ⁿ)
1 2 6 7 10 14 15 16 21 22 28 29 30 31 32 33 34	4 5 4 4 5 6 6 6 6 7 7 7 8 6 6 7 7	0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4	9.950120(0) 9.801984(0) 9.559971(0) 9.231160(0) 8.824965(0) 8.352683(0) 7.827014(0) 7.261469(0) 6.669735(0) 6.065274(0) 5.460716(0) 4.867497(0) 4.295553(0) 3.753093(0) 3.246508(0) 2.780359(0) 2.357449(0) 1.978977(0) 1.644736(0) 1.353346(0) 1.102500(0) 8.892118(-1) 7.100501(-1) 5.613449(-1) 4.393672(-1) 3.404729(-1) 2.612128(-1) 1.984100(-1) 1.492071(-1) 1.110894(-1) 8.188665(-2) 5.975999(-2) 4.317823(-2) 3.088704(-2)	9.95012479(0) 9.80198673(0) 9.55997478(0) 9.23116344(0) 8.82496903(0) 8.35270199(0) 7.82704545(0) 7.26149030(0) 6.66976789(0) 6.06530660(0) 5.46074412(0) 4.86752228(0) 4.29557318(0) 3.75311111(0) 3.24652467(0) 2.78037290(0) 2.35746057(0) 1.97898674(0) 1.64474464(0) 1.35335283(0) 1.10250492(0) 8.89216081(-1) 7.10053615(-1) 5.61347500(-1) 4.39369336(-1) 3.40474421(-1) 2.61214065(-1) 1.98410974(-1) 1.49207819(-1) 1.11089965(-1) 8.18869738(-2) 5.97602198(-2) 4.31784069(-2) 3.08871441(-2)		35 : : : : :50 51 52 53 54 55 56 57 58 59 60 61 : : :88 89 90 91	7 : : : : :10 10 11 11 11 11 11 11 11 11 11 11 : : :11 10 11 11	3.5 : : : : :5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0 6.1 : : :8.8 8.9 9.0 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10.	2.187483(-02) : : : : :3.726643(-05) 2.249050(-05) 1.343809(-05) 7.949373(-06) 4.655704(-06) 2.699572(-06) 1.549749(-06) 8.808157(-07) 4.956393(-07) 2.761235(-07) 1.522994(-07) 8.316679(-08) : : :1.527919(-16) 6.306696(-17) 2.576850(-17) 1.042381(-17) 4.175947(-18) 1.660551(-18) 6.670604(-19) 3.025806(-19) 2.349915(-19) 4.096942(-19) 1.019628(-18) 2.731899(-18) 7.452993(-18)	2.18749112(-02) : : : : :3.72665317(-05) 2.24905706(-05) 1.34381028(-05) 7.94938558(-06) 4.65571332(-06) 2.69957850(-06) 1.54975396(-06) 8.80816483(-07) 4.95639984(-07) 2.76124090(-07) 1.52299797(-07) 8.31670729(-08) : : :1.52797740(-16) 6.30615778(-17) 2.57675711(-17) 1.04240256(-17) 4.17501738(-18) 1.65551973(-18) 6.49931301(-19) 2.52616378(-19) 9.72094942(-20) 3.70353883(-20) 1.39694133(-20) 5.21670711(-21) 1.92874985(-21)

. . When IB is 91, the solutions are regarded as being convergent at all the points after t=9.1, because there are not eleven equidistant points in the range 9.1≤t≤10 if IB is increased to 92. The solutions after t=9.5 are very different from the exact solutions but these are useful to make the solution convergent at t=9.1. The case occurs always when the IB value is increased only by one as shown in Table 1.6. When the solution is convergent at t=5.2, the solutions after t=5.2 are not convergent and the solution at t=6 is as inaccurate as the result in R.-K. 3 of Table 1.1. It is said that by this reason Milne's method uses the predictor introduced by the backward interpolation formula and that it is able to predict the solution accurately at the last point. However, the idea is a mistake because Milne's predictor is of the same as Eq. (1.14) introduced by the forward interpolation formula. It is also said that the predictor-corrector methods using the polynomial of higher order are unstable. However, the cause is that the methods predict and correct the solution at the last point.

Table 1.6 . .From 0 to 6 by h=0.1.
IB	E	t	y (10⁻ⁿ)	True (10⁻ⁿ)
: :45 46 47 48 49 50 51 52	: :9 9 9 10 10 10 10 10	: :4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0	: :4.006519(-4) 2.541928(-4) 1.596674(-4) 9.929478(-5) 6.113552(-5) 3.726643(-5) 2.249050(-5) 1.343809(-5) 7.949374(-6) 4.655708(-6) 2.699582(-6) 1.549770(-6) 8.808545(-7) 4.957133(-7) 2.762749(-7) 1.526310(-7)	: :4.00652974(-4) 2.54193577(-4) 1.59667624(-4) 9.92949522(-5) 6.11356511(-5) 3.72665317(-5) 2.24905706(-5) 1.34381028(-5) 7.94938558(-5) 4.65571332(-6) 2.69957850(-6) 1.54975396(-6) 8.80816483(-7) 4.95639984(-7) 2.76124090(-7) 1.52299797(-7)

. . The numerical methods to solve the differential equations with ill condition must satisfy the following conditions. If they are of using the solution by interpolation, the equidistant interval must be variable and the integral matrix may be of variable order or fixed order. They must decide at the last point whether the solution is convergent or not. If they are of using the solution by extrapolation, the integral matrix is of variable order. They must decide at the next point of the initial value whether the solution is convergent or not. If the solutions are convergent at the successive k points, the next step may advance the initial point to the k-th point.

1.5 Taylor's expansion of a vectorial function.
. . In order to solve the differential equation by the author's operational calculus using vector computer, the vector-valued function F(Y) of the differential equation y' =f(y) must be calculated by use of the vector Y immediately. In case of the usual computer, the vector values may be obtained by calculating the successive differences of the values f_k for the solutions y_k transformed from the vector Y at equidistant points. In the case, if the calculation causes the successive differences the cancellation of many significant figures, the solution may not get sufficient accuracy. In order to avoid it, the vector-valued function should be calculated by use of the vector Y immediately. Accordingly, the function F(Y) must be the vectorial function mentioned in Section 1.2. If the function is an irrational function it is calculated by use of Taylor's expansion usually.
. . Carrying out Taylor's expansion of the function f(t) at every equidistant point t_k=t₀+kh for k≥0, and calculating the difference of the initial value f(t_k) from the value f(t_k+h) at t=t_k+h of the every Taylor's expansion,

Calculating the successive differences of Δf₀, Δf₁, Δf₂, Δf₃, ………,
Eq1_29c

Hence these Taylor's expansions are expressed in vector formula as in Eq. (1.30).

. . [Theorem 401]. .When f(t)∈C^p+1[a, b], Taylor's expansion of the vector ΔF which denotes the function Δf(t) is expressed in,

_____(1.30) [Proof]. .Δ^kf(t) satisfies Δ^kf(t)∈C^p+1[a, b] because of f(t)∈C^p+1[a, b] and,
Eq1_30a

Taylor's expansion of Δ^kf(t) at t₀ gives the following at t=t₀+h.
Eq1_30b

Because Δ^k+1f(t₀) and Δ^kf^(m)(t₀) for k≥0 and 1≤m≤p are the elements of the vectors F and F^(m), the equations (A) for k≥0 are expressed in Eq. (1.30) with the vectors.________[Q.E.D.]

. . Accordingly, if the function f(y) is an irrational function of y, the differential equation Y' =F is equivalent to the differential equation of the (p+1)th order by substituting it into F of Eq. (1.30). However, if the successive derivatives of f(y(t)) is independent of y, the differential equation can be solved as the differential equation of the first order. The example is the differential equation in Eq. (1.18) of Section 1.3.

. .(1.31) This is equivalent to solving the differential equation whose right side is Taylor's expansion of f(t) as follows.
Eq1_31a

. . The function f(t) is expressed in the vectorial function F(T) by replacing the denotation f and t with the respective vectors F and T as mentioned in Chapter 3 Section 1.7. Accordingly, Taylor's expansion of f(t) is expressed in the vectorial function as follows.
Eq1_32

____(1.32). . If the variable t takes the equidistant points t₀, t₁, t₂, ………, t_m, the function f(t) is expressed in the vector whose elements are the initial value f₀=f(t₀) and the successive differences and it is denoted by the vector F₀. When the initial value is the value f(t) for any point t, the function is denoted by the vector F. If the points t₀, t₁, t₂, ………, t_m are not always equidistant, the variable t is a function of another variable u denoted by t(u). If the variable u takes the equidistant points u₀, u₁, u₂, ………, u_m, the function t(u) is expressed in the vector denoted by T₀. When the initial value is the value t(u) for any point u, the function t(u) is denoted by the vector T. In the case, the function f(t) is expressed in the vectorial function F(T) because f(t) is the function of the function t(u). These are the meaning of every denotation in Taylor's expansion of the vectorial function F(T) expressed in Eq. (1.32).
. . When the vector T₀ denotes the function t(u) and the vector T denotes the function t(u)+h, the vector T consists of the successive differences of t₀+h, t₁+h, t₂+h, ………, t_m+h. Hence,
Eq1_32a

Accordingly, Eq. (1.32) becomes Eq. (1.30), where all the subscript zeros are omitted and the points t₀, t₁, t₂, ………, t_m are equidistant by the interval h.
. . When the vector T denotes the function t(u) from u₀ to u_m and the vector T₀ denotes the constant function t=t(u₀)=t₀,
Eq1_32b

____(1.33) The coefficients f₀^(k) can be obtained by usual differentiation of the function f(t). Hence, the value of the vectorial function F(T) can be obtained by Eq. (1.33) for given vector T.
. . When the function t(u) is the function such that t(u)=t₀+hu, the value of t(u) for the equidistant points u=0, 1, 2, 3, ……… are also equidistant and the function t(u) is expressed in the vector T=(t₀, h, 0, 0, ………). Hence,
Eq1_33a

Substituting these vectors into Eq. (1.33) and omitting the first element f₀ of the vectorial function F(T) because it is the value given and all the first elements of these vectors are zero, Eq. (1.33) is expressed in the matrix equation as follows.
Eq1_34

__(1.34). . The equation also can be obtained by use of the values at the equidistant points t₀, t₁, t₂, ……… of Taylor's expansion of the function f(t) at t=t₀.
Eq1_34a

The successive differences of these equations give zeros to all the first terms of the right sides and give the following coefficient to the term of h^mf₀^(m) of the difference Δⁿf₀.
Eq1_35

___(1.35)
where, _____n, m=1, 2, 3, 4, ………. In case of n=1, the coefficients a(1, m) become the first row of the matrix of Eq. (1.34). By Theorem 4.21 of Chapter 3 Section 4, all the coefficients a(n, m) are one in case of m=n and zero in case of m<n. These are equal to the diagonal elements and the under elements of them. In case of m>n, the coefficients a(n, m) are obtained by the following theorem and equal to the upper elements of the diagonal.

. . [Theorem 402]. .In case of m≥n+1,

______(1.36) [Proof]. .Because of m−1≥n,
Eq1_36a

[Q.E.D.]

. . In case of the differential equation in Eq. (1.31), all the values f₀", f₀⁽³⁾, f₀⁽⁴⁾, ……… are zero, so the matrix equation Eq. (1.34) becomes,
Eq1_36b

This is of the same as the integrand F in Eq. (1.31).
. . In general, if the differential equation is differentiable, the initial values y₀^(k) for the successive integer k can be obtained and the solution y(t) can be obtained by Taylor's expansion without integration. The numerical solutions can also be obtained by the matrix equation in Eq. (1.34), where the denotation F and f must be replaced with Y and y respectively. Accordingly, it is equivalent to the convergence of Taylor's expansion that there exists the numerical solution of the differential equation.
. . When we solve the differential equation of the first order expressed in y' =f(y) by use of the iterative integration, if the function f(y) is irrational, it may be expressed in Taylor's expansion which is the function of the independent variable y. The vectorial function F(Y) is expressed in Eq. (1.37) by replacing the denotation T and t of Eq. (1.33) with Y and y respectively.
Eq1_37

____(1.37) The values y₀, y₁, y₂, ……… are not always equidistant. Accordingly, the vector Y−y₀ is,
Eq1_37a

Expanding the vector into the matrix,
Eq1_37b

Therefore, the value of the vectorial function F(Y) can be calculated by use of the vector Y.
. . When the functions y(t) and f(y) are,
Eq1_38

______(1.38) the value of the vectorial function F(Y) can be calculated with the vector Y as follows.
Eq1_38a

Differentiating f(y) by the variable y and substituting y₀ into the derivatives,
Eq1_39

__(1.39). . The result is also obtained by calculating the successive differences of the values of the function f(y(t)) at the equidistant points t_k=kh for k≥0. Differentiating f(y(t)) by t,
Eq1_39a

Accordingly, the Taylor's expansion of the function f(y(t)) is,
Eq1_40

____(1.40) The values at t=0, h, 2h, 3h, ……… are,
Eq1_40a

The successive differences are of the same as the elements of the vector F(Y) in Eq. (1.39).

1.6 Newton-Raphson's method.
. . Denoting the square root of a number y by x, it is given by the point at which the locus of the function f(x)=x²−y intersects the x-axis. In the same way, denoting the square root of a vector Y by the vector X, it is given by the point at which the locus of the vectorial function F(X)=X²−Y intersects the X-axis on the plane denoting all the ordered pairs of vectors (X, F). The point is obtained by Newton-Raphson's method.

_____(1.41)
Eq1_42

___(1.42). . When the function y(t) is the quadratic function in Eq. (1.38), the iterative method may be started on the assumption x₀(t)=1, that is, the square root of y₀. These functions are expressed in the vectors as follows.
Eq1_42a

Denoting the quotient Y/X by W, the quotient W is obtained by separating the diagonal elements from the matrix expansion X of the vector X in the equation Y=WX =XW by Theorem 1.6 of Chapter 3 Section 1.6 as expressed in Eq. (1.43). The first element w₀ is obtained first as y₀/x₀ independently of the vector W on the right side. The second element Δw₀ is obtained by substituting the obtained w₀ into the right side. The third element Δ²w₀ is obtained by substituting the obtained Δw₀ into the right side and so on. Accordingly, the quotient W has an infinity of elements generally. However, it may be truncated into the vector which has the successive differences up to the order needed.
Eq1_43

____(1.43) When the vector X is the vector 1, all the elements of the matrix of Eq. (1.43) are zero and the quotient W is equal to the vector Y. Hence,
Eq1_43a

The second iteration becomes,
Eq1_43b

. . The result X₂ is equal to Taylor's expansion in Eq. (1.40) up to the term of h⁶. Newton-Raphson's method is a better method to obtain the square root of a vector than Taylor's expansion, because it does not need the successive derivatives of the square root of the function y(t).

1.7 The formulas for differentiation and integration.
. . The author's differential and integral operators may be used as the symbolical operators to solve the differential equations. In the case, it must be cautioned that the multiplication of the operators and vectors does not always hold the same commutative and associative laws as the multiplication of numbers as mentioned in Chapter 2 Section 2.2. The differential equation in Eq. (1.18) of Section 1.3 may be solved by no use of the square root of the function y, squaring the both sides.

________(1.44) The differential equation must not be changed as follows.
Eq1_44a

The change is not admitted by the associative law of the multiplication of the operator and the product of two vectors 1•Y' , because it is equivalent to the change of the integral operation as,
Eq1_44b

The differential equation in Eq. (1.44) can be solved correctly by use of the formulas which are equivalent to the formulas for differentiation and integration of usual functions.

. . [Theorem 403]

________(1.45) [Proof]. .The derivative of yⁿ(t) is,
Eq1_45a

Denoting the values of the derivative at equidistant points t_k=t₀+kh by (yⁿ)' _k, the both sides are expressed at t=t₀+qh for the real number q in Newton's interpolation formula as,
Eq1_45b

At all the equidistant points,
Eq1_45c

Therefore,
Eq1_45d

By the definition of the product of two vectors,
Eq1_45e

____[Q.E.D.]. . Differentiating the differential equation in Eq. (1.44) by Theorem 403, it is solved as,
Eq1_45f

Carrying out the above integral by use of the vector expressed in the elements,
Eq1_45g

Theorem 403 is held also when the n value is not integer. Hence the differential equation in Eq. (1.18) may be solved as follows.
Eq1_45h

. . The differential equation y' =y⁻¹, y(t₀)=y₀ also is solved as follows.
Eq1_45i

. . When the integrand is the derivative f(t) of a function g(t) with the initial value g₀=g(t₀), denoting these functions in the vectors F and G,
Eq1_46

____(1.46) Accordingly, the operator hΔ^TD⁻¹ is equivalent to the integral operator below.
Eq1_47

________(1.47) The integral operator h(ΔD)⁻¹ adds the initial value g₀1 to the result in Eq. (1.46) and obtains the solution of the differential equation g' =f(t). If the initial value is unknown it may add an arbitrary constant denoted by c1, supposing g₀=0 and it obtains the indefinite integral.
. . The indefinite integral does not mean that there are many solutions of the differential equation but that the initial value is unknown. Hence the integral operation is unique. It is said that the differential operation obtains the one derivative g' (t) from all the functions with any constant c, that is, g(t)+c. However, the author says that the differential operation obtains the pair of the initial value and the derivative, that is, (g(0), g' (t)) and that the two operations are reversible.

. . [Theorem 404]

______(1.48) [Proof]. .The derivative of the product of two functions y(t)z(t) is (yz)' =y' z+yz' . Expressing the equation in Newton's interpolation formula, it is expressed in the vector equation in Eq. (1.48) by the same way as the proof of Theorem 403._________________________________[Q.E.D.]

. . [Theorem 405]

________(1.49) [Proof]. .It is evident by integrating the both sides of Eq. (1.48).__________[Q.E.D.]

. . The differential equation y' =y⁻¹, y(t₀)=y₀ may be solved by Theorem 405 as follows.

This solution Y is the expression in the function of the vector T, that is, the vectorial function. When the vector T is (t₀, h, 0, ……) the vector Y denotes the function y(t) which is the solution of the differential equation y' =y⁻¹, y(t₀)=y₀. Even if the vector T is the other values denoted by (t₀, Δt₀, Δ²t₀, Δ³t₀, ……) the vector Y satisfies the differential equation Y' =Y⁻¹ and the formulas for differentiation and integration. Hence the formulas for differentiation and integration are held too in case that the two vectors Y and Z are vectorial functions Y(T) and Z(T). In the case, the differential operator is denoted by d/dT and the inverse operator is denoted by,
Eq1_49b

When the vectorial function Y(T) gives the vector for a vector T, the author may denote these operators simply by,
Eq1_49c

When the vector T is (t₀, h, 0, ……) the operators are denoted by h⁻¹ΔD and h(ΔD)⁻¹ respectively. In the case that the vector T is an arbitrary value expressed in (t₀, Δt₀, Δ²t₀, Δ³t₀, ……), it is able to use these operator only as the symbols of operation. However, if the vector T is transformed into the vectorial function of the vector U whose value is (u₀, h, 0, ……) denoting equidistant interval, it is able to use them as the matrix.

. . [Theorem 406]. .When Y=e^T, the derivative is Y' =e^T.
[Proof]. .It is evident by y=e^t and y' =e^t.__________[Q.E.D.]

. . The vectorial function e^T is expressed in Taylor's expansion as,

________(1.50) It is evident that the differentiation of the expansion by d/dT gives e^T. When the vector T is (0, h, 0, ……), by the equations before Eq. (1.34),
Eq1_50a

Accordingly, the value of the vectorial function F(T)=e^T for the T value is given as,
Eq1_51

____(1.51) Premultiplying the vector F by the differential operator h⁻¹ΔD,
Eq1_51a

The other formulas for differentiation and integration are all held as those of the vectorial functions and the vectors.

Return to CONTENTS