Note that the factorization below
suggests a the simpler equation,
i.e. translating $(-2,-3)$ to the origin
$$ y'=(x+2)(y+3)-6 $$
and the initial condition $(0,-1)$ to $(2,2)$.
This, however, leaves $(y+3)'=y'$ unchanged:
$$
\eqalign{
y' &= xy - 6 \\
y' - xy &= - 6 \\
(y' - xy) \, e^{-x^2/2} &= - 6\,e^{-x^2/2} \\
(y \, e^{-x^2/2})' &= - 6\,e^{-x^2/2} \\
y \, e^{-x^2/2} &= - 6 \int e^{-x^2/2} dx \\
y \, e^{-x^2/2} &= c-6 \sqrt{\frac\pi2}\,
\text{erf}\frac{x}{\sqrt{2}}
\qquad \implies \quad c=y_0 \\
y &= e^{x^2/2} \left( y_0 - 3 \sqrt{2\pi}\,
\text{erf}\frac{x}{\sqrt{2}}\right) .
}
$$
Translating back, we get
$$
y = e^{(x+2)^2/2} \left( c - 3 \sqrt{2\pi}\,
\text{erf}\frac{x+2}{\sqrt{2}}\right)-3
$$
for the same constant $c$ which, however,
is no longer equal to $y_0$, but rather to
$$
\eqalign{
c &= (y_0+3) e^{-\frac12(x_0+2)^2}
+ 3\sqrt{2\pi}\,\text{erf}\frac{x_0+2}{\sqrt{2}}
\\&=2e^{-2}+3\sqrt{2\pi}\,\text{erf}\sqrt{2}
\approx 7.44839864640532
}
$$
which would give
$$
\eqalign{
y(0.1)
&= e^{1/200} \left( c - 3 \sqrt{2\pi}\,
\text{erf}\frac{2.1}{\sqrt{2}}\right)-3
\\&\approx -1.21143171766941
\,.
}
$$
If we had just let $u=x+2$ and $v=y+3$ to begin with,
we could have rewritten our differential equation as
$v'=uv-6$
which is easier to differentiate by hand anyway:
$$
\eqalign{
\\v^{(1)}&=\left(u\right)\cdot&v+\left(-6\right)
\\v^{(2)}&=\left(u^{2}+1\right)\cdot&v+\left(-6 \, u\right)
\\v^{(3)}&=\left(u^{3}+3\,u\right)\cdot&v+\left(-6\,u^{2}-12\right)
\\v^{(4)}&=\left(u^{4}+6\,u^{2}+3\right)\cdot&v+\left(-6\,u^{3}-30\,u\right)
\\v^{(5)}&=\left(u^{5}+10\,u^{3}+15\,u\right)\cdot&v
+\left(-6\,u^{4}-54\,u^{2}-48\right)
\\v^{(6)}&=\left(u^{6}+15\,u^{4}+45\,u^{2}+15\right)\cdot&v
+\left(-6\,u^{5}-84\,u^{3}-198\,u\right)
\\v^{(7)}&=\left(u^{7}+21\,u^{5}+105\,u^{3}+105\,u\right)\cdot&v
+\left(-6\,u^{6}-120\,u^{4}-522\,u^{2}-288\right)
\\v^{(8)}&=\left(u^{8}+28\,u^{6}+210\,u^{4}+420\,u^{2}+105\right)\cdot&v
+\left(-6\,u^{7}-162\,u^{5}-1110\,u^{3}-1674\,u\right)
\\v^{(9)}&=\left(u^{9}+36\,u^{7}+378\,u^{5}+1260\,u^{3}+945\,u\right)\cdot&v
+\left(-6\,u^{8}-210\,u^{6}-2070\,u^{4}-5850\,u^{2}-2304\right)
}
$$
especially if one notices the recursion
$$
\left.\matrix{
v^{(n)}=a_nv+b_n\\\\
a_1=u,~b_1=-6}\right\}
\quad\implies\quad
\left\{\matrix{
a_{n+1}=a_n'+u\,a_n\\\\
b_{n+1}=b_n'-6\,a_n\,,}\right.
$$
which I managed to calculate in sage:
var('a,b,u,v,z')
a(u) = u
b(u) = -6
for n in range(1,10):
# print n, a, b #Note: a & b are functions of u
# print n, a(u).simplify_full(), b(u).simplify_full()
print '\\\\v^{(%d)}&=\\left(%s\\right)\cdot v+\\left(%s\\right)' \
% (n, latex(a(u).simplify_full()), latex(b(u).simplify_full()))
# print 'v^{(%d)}=%d' % (n, 2*a(2)+b(2))
z(u) = a
a(u) = (diff(a,u) + u*z).simplify_full()
b(u) = (diff(b,u) - 6*z).simplify_full()
Using the initial condition $(u,v)=(2,2)$
corresponding to $(x,y)=(0,-1)$, we get
$$
v^{(1)}=-2\\
v^{(2)}=-2\\
v^{(3)}=-8\\
v^{(4)}=-22\\
v^{(5)}=-76\\
v^{(6)}=-262\\
v^{(7)}=-980\\
v^{(8)}=-3794\\
v^{(9)}=-15428
$$
so that, factoring out the $(-1)$,
I computed the following terms and partial sums:
t = [1, 2, 2, 8, 22, 76, 262, 980, 3794, 15428]
x = 0.1
y = 0
for n in range(10):
tn = t[n]*x^n/factorial(n)
y = y + tn
print '%-4d %-14g %-16.12g' % (n, tn, y)
0 1 1
1 0.2 1.2
2 0.01 1.21
3 0.00133333 1.21133333333
4 9.16667e-05 1.211425
5 6.33333e-06 1.21143133333
6 3.63889e-07 1.21143169722
7 1.94444e-08 1.21143171667
8 9.40972e-10 1.21143171761
9 4.25154e-11 1.21143171765
So the analytic and Taylor series solutions agree.
To be complete, it is worth noting that there are many variants of Runge-Kutta (RK), including some new ones. These include, but are not limited to:
- Generalized RK, RK1, RK2, RK3, RK4, Simple RK, Euler-Cauchy, Optimal Method, RK5, New Variant RK, RK-Merson (RKM), RK-Fehlberg (RKF), Cash-Karp-RK (CKRK), Implicit RK (IRK) ...
For this problem, I will assume you want to typical fourth order RK (RK4).
Initialization
- $h = \dfrac{b-a}{N} = .1 = \dfrac{2-1}{N} \rightarrow N = 10$
- $x_0 = a = 1, y_0 = y(a) = y(1) = \alpha = 0 \rightarrow a = 1, \alpha = 0$
Iteration using RK4
- $y'(x) = f(x, y) = \dfrac{3 x -2 y}{x}, y(1) = 0$
- $k_1 = h f(x_i, y_i)$
- $k_2 = h f(x_i + h/2, y_i + k_1/2)$
- $k_3 = h f(x_i + h/2, y_i + k_2/2)$
- $k_4 = h f(x_i + h, y_i + k_3)$
- $y_{i + 1} = y_i + \dfrac{1}{6}\left( k_1+2k_2+2k_3+k_4 \right)$
- $x_{i + 1} = x_i + h = 1 + 0.1 i$
Results
- For $i= 1$, we have:
- $x_0 = 1, y_0 = 0$
- $k_1 = h f(x_i, y_i) = \dfrac{1}{10} \left(\dfrac{3 x_{i} - 2 y_{i}}{x_{i}}\right) = \dfrac{1}{10}\left(\dfrac{3 x(0) - 2y(0)}{x(0)}\right) = \dfrac{1}{10}\left(\dfrac{3 (1) - 2(0)}{1}\right) = 0.3$
- $k_2 = h f(x_i + h/2, y_i + k_1/2) = \dfrac{1}{10} \left(\dfrac{3 (x_{i}+(1/20)) - 2(y_{i} + (0.3/2))}{x_{i}+(1/20)}\right) = \dfrac{1}{10} \left(\dfrac{3 (1+(1/20)) - 2(0 + (0.3/2))}{1+(1/20)}\right) = 0.271428571428571428571428571428571428571428571428571428571428$
- $k_3 = h f(x_i + h/2, y_i + k_2/2) = \dfrac{1}{10} \left(\dfrac{3 (x_{i}+(1/20)) - 2(y_{i} + (0.27142857/2))}{x_{i}+(1/20)}\right) = \dfrac{1}{10} \left(\dfrac{3 (1+(1/20)) - 2(0 + (0.27142857/2))}{1+(1/20)}\right) = 0.274149659863945578231292517006802721088435374149659863945578$
- $k_4 = h f(x_i + h, y_i + k_3) = \dfrac{1}{10} \left(\dfrac{3 (x_{i}+(1/10)) - 2(y_{i} + 0.2741496598)}{x_{i}+(1/10)}\right) = \dfrac{1}{10} \left(\dfrac{3 (1+(1/10)) - 2(0 + 0.2741496598)}{1+(1/10)}\right)= 0.250154607297464440321583178726035868893011750154607297464440$
- $y_{i + 1} = y_i + \dfrac{1}{6}\left( k_1+2k_2+2k_3+k_4 \right) = 0.273551844980416408987837559266130694702123273551844980090476 \approx 0.273552$
- $x_{i + 1} = x_i + h = 1 + 0.1 i = 1.1$
Continuing with this itearation, we generate the following table.
$~~~~~\text{Step} ~~|~~~~ x ~~~|~~ y $
- $~~00 ~~|~~ 1.0 ~~|~~ 0. 00000$
- $~~01 ~~|~~ 1.1 ~~|~~ 0.273552$
- $~~02 ~~|~~ 1.2 ~~|~~ 0.505553$
- $~~03 ~~|~~ 1.3 ~~|~~ 0.708281$
- $~~04 ~~|~~ 1.4 ~~|~~ 0.889793$
- $~~05 ~~|~~ 1.5 ~~|~~ 1.05555$
- $~~06 ~~|~~ 1.6 ~~|~~ 1.20937$
- $~~07 ~~|~~ 1.7 ~~|~~ 1.35398$
- $~~08 ~~|~~ 1.8 ~~|~~ 1.49136$
- $~~09 ~~|~~ 1.9 ~~|~~ 1.62299$
- $~~10 ~~|~~ 2.0 ~~|~~ 1.75$
The exact solution is given by:
$$y(x) = \dfrac{x^3-1}{x^2}$$
At $x=2$, we have: $y(2) = \dfrac{7}{4} = 1.75$
Compare that to Runge-Kutta's-4 method, which has $1.75$.
Best Answer
So it turns out, as happens, that it was a mathematical error in the implementation that yielded a very small absolute error, related to the presence and absence of the time step multiplier on several terms.
The first error was to miss multiplying the small state deltas in the second, third and fourth tangent calculations for the second state. More importantly I had combined two varying definitions of the method, resulting in one too many multiplications by the time step when adding the slope delta to the current state values.
The warning here I suppose is to ensure a consistent set of equations is sourced while developing the method. I had drawn from several varying sources.