2

I'm using a Runge-Kutta second-order solver (more specifically, this with $\alpha = 1$) and I'm trying to approximate the solution of the nonlinear system

$$\begin{cases} T'(r)=S(r) \\ S'(r)=-\frac{S(r)}{r}-\beta e^{T(r)} \end{cases}$$

with initial values $T(1)$ and $S(1)$. Moreover, I'm interested in the solution in the interval $[h, 1]$, where $h$ is the step size used in the solver. Notice that to go back in time I'm using a negative step $h_-=-h$ to run the algorithm. The solutions look correct and I'm computing the convergence order with the formula

$$\log_2\left(\frac{\text{mean}(|T_h - T_{h/2}|)}{\text{mean}(|T_{h/2} - T_{h/4}|)}\right)$$

where $h$ is the step size used to compute the solution $T$ (I'm only interested in finding $T$, not $S$, for this specific problem). This expression gives values of about 2 for the first 4 values of $h$ starting from $0.05$, although they are dropping to $\approx 1.95$ after 4 iterations. After $15$ iterations, it seems to stabilize to $\approx 1$.

I'm fairly certain that this is not due to the implementation of the method since I tried the same code with other non-singular at $r=0$ systems and, with those systems, obtained values approaching the correct (theoretical) value of $2$.

Is there anything more to this than roundoff errors in the division $\frac{S}{r}$ in the $S'$ term? When $h$ is very small then $r$ can be very small near the origin. Can this behaviour be somewhat diminished while still using this R-K method to solve the system of equations?

Edit: one more interesting thing about this problem is that if we only compute the convergence order using values of $T_h$ far from the singularity, we obtain the correct convergence order of $\approx 2$.

Edit 2 (more information): I didn't want to give out too many details about this problem, but here goes the full thing. I'm interested in solving the BVP

$$T''(r) + \frac{1}{r}T'(r) + \beta e^{T(r)} = 0$$

with boundary conditions $T(1)=0$ and $T'(0)=0$. Originally I used the Newton method together with the power series approximation $T(h) - T(0) = -\frac{\beta e^{T(0)}}{4}h^2$ to avoid the singularity and this gave me a convergence order of $2$. Now I'm interested in solving this with a shooting method using the trapezoidal method which I mentioned in the first paragraph of the post.

The difficulty here is dealing with the singularity, hence I solved the problem in $[h, 1]$ using a negative step as I described above. Now I'll describe my current algorithm:

  1. let $s_1$, $s_2$ be such that $S(0) > 0$ for $S(1) = s_1$ and $S(0) < 0$ for $S(1) = s_2$.
  2. compute $T(0)$ using the power series approximation for $S(1) = \frac12 (s_1 + s_2)$. That is, find $T(0)$ using Newton's method with a specified tolerance, I used $10^{-12}$. This converges very fast since the starting value I'm giving the method is $T(h)$.
  3. approximate $T'(0)$ using a second order forward finite difference (that is, depending on $T(0)$, $T(h)$ and $T(2h)$). I call this a "compatibility condition", since the power series forces the first order finite difference to be very close to $0$.
  4. depending on the sign of the approximation obtained, rerun the program with $\frac12 (s_1 + s_2)$ as one of either $s_1$ or $s_2$. Essentially, we're using the bisection method over $[s_1, s_2]$ to find a value of $s$ that gives good values for $T'(0)$. Stop running if $s_2 - s_1$ is below a specified tolerance (e.g. $10^{-12}$), and classify the average of these two as the correct approximation for $s$.

I've worked on this for some time after originally posting this question, and this method seems to be giving good (and fast) convergence. The values I'm getting while using the formula given above in the post are all around $2$.

There might exist some theoretical problems with this method, but it seems to work well in practice.

Originally, I was experiencing problems with the convergence order running the program with a fixed $s$ right from the first value of $h$ for which I computed an approximation. I figured this can be partly explained by the paper I posted in one of my comments.

user401936
  • 1,071
  • Thank you for this problem. I do not have time to replicate the calculations right now. I am not certain that taking the mean of the error will give meaningful results. I am, however, certain that the fractions $F_h = \frac{T_{2h}(t) - T_{4h}(t)}{T_h(t) - T_{2h}(t) }$ will (eventually) converge monotonically towards $2^p$ for fixed $t$ as you reduce $h$ to zero. Here $p=2$ for Heyn's method. Now, this is in exact arithmetic only. The computed values will follow the same pattern until subtractive cancellations sets in. – Carl Christian Jan 13 '20 at 08:36
  • I now have a working example with demonstrates that the order of accuracy is $p=2$ for all sample points except the singularity $r=0$. I need your exact value of beta, as well as your initial conditions to make further progress. – Carl Christian Jan 13 '20 at 17:10
  • Hello Carl, that's interesting. I'm running a value of $\beta = 0.3$ and initial conditions $T(1) = 0$ and $S(0) = -0.15609052263065767$. More interestingly, I found a paper that seems to analyze the behaviour I'm observing: https://www.jstor.org/stable/2007795?seq=1

    Notice that example 1 in section 4 is this exact problem! Now I'm even more confused as to how you managed to get the correct convergence orders. The first 10 values I'm getting are [2.00174681, 1.9968536, 1.98806133, 1.97286816, 1.94525955, 1.89537815, 1.81062822, 1.68245844, 1.51943556, 1.35336385].

     – user401936 Jan 13 '20 at 17:53
  • Also, other starting conditions only give me a convergence order of 1. I should add that the value $S(0)$ is being obtained by an approximation using a shooting method -- it seems that the convergence is very dependent on picking a good value of $S(0)$. – user401936 Jan 13 '20 at 17:59
  • Note: your original question uses S and T at r=1, i.e., the same point. – Carl Christian Jan 13 '20 at 18:26
  • My mistake, $S(0)$ in the above comments is in fact $S(1)$. – user401936 Jan 13 '20 at 18:27
  • I am pressed for time right now, but look forward to returning to this question. I believe that it is crucial that you also treat the value of $T(0)$ as a function of $h$. The integration of each shot must be done with a method which has order at least 2. Moreover, you must use a very small tolerance your shooting method. – Carl Christian Jan 14 '20 at 20:25
  • I finished my other work. There is still some lingering doubt about what you have and want. It seems to me you are given $S(1)$ and must find $T(1)$ such that $S(0)=0$. Is this correct? Please do not confirm this in a comment, but update your question instead. Include all the specific values that you are interested in. I am having trouble getting the shooting method to converge. – Carl Christian Jan 20 '20 at 15:33
  • I added a lot of details to the post. Originally I missed a lot of important information since I didn't want to give too much away, but more details are necessary than I had anticipated. I'm not too sure about posting my Python code online, but since you are very interested in this problem and want to see the current code I could find a way of giving it to you.

    If my current edit is still not very clear, I am given $T(1)$ and must find $S(1)$ s.t. $S(0) = 0;(=T'(0))$.

     – user401936 Jan 21 '20 at 16:45
  • No. That's a mistake in my post, which I'll edit. – user401936 Jan 22 '20 at 03:54

1 Answers1

1

Note that applying some limit formula to the middle term for $r\to 0$ results in $2T''(0)+βe^{T(0)}=0$, so that the Taylor formula of degree $2$ for $T(h)$ is $$ T(h)\approx T(0)-\fracβ4e^{T(0)}h^2. $$ Simple forward integration from $r=h$ using $S(h)=T'(h)=-\fracβ2e^{T(0)}h$ and finding the root $T(0)$ for $T(1)=0$ with the secant method gives the table $$\small \begin{array}{ll|llll} N & h=1/N & T_h & T_h-T_{h/2} & \frac{T_h-T_{h/2}}{T_{2h}-T_h} & -\log_2(...) \\ \hline 10 & 0.1 & 0.079645659009 & 2.776495e-05 & 0.240477 & 2.056027 \\ 20 & 0.05 & 0.079617894057 & 6.625231e-06 & 0.238618 & 2.067222 \\ 40 & 0.025 & 0.079611268826 & 1.602851e-06 & 0.241931 & 2.047330 \\ 80 & 0.0125 & 0.079609665974 & 3.929955e-07 & 0.245185 & 2.028056 \\ 160 & 0.00625 & 0.079609272979 & 9.720956e-08 & 0.247355 & 2.015343 \\ 320 & 0.003125 & 0.079609175769 & 2.416719e-08 & 0.248609 & 2.008048 \\ 640 & 0.001563 & 0.079609151602 & 6.024530e-09 & 0.249285 & 2.004129 \\ 1280 & 0.0007813 & 0.079609145578 & 1.503949e-09 & 0.249638 & 2.002093 \\ 2560 & 0.0003906 & 0.079609144074 & 3.757122e-10 & 0.249817 & 2.001055 \\ \end{array} $$ where none of the reported difficulties in the convergence rate are visible, the empirical order is stably $2$.

Computing the next term gives $$ T(h)\approx T_4(h)=T(0)-\fracβ4e^{T(0)}h^2+\frac{β^2}{64}e^{2T(0)}h^4. $$

Now using the classical RK4 method from $r=h$, using the approximate values $T_4(h)$ and $T_4'(h)$ as initial values, and applying the secant method to find the value of $T(0)$ so that $T(1)=0$ gives the table \begin{array}{ll|ll} N&h=1/N&T(0)&\text{error estimate}\\ \hline 5 & 0.2 & 0.07960889449312651 \\ 10 & 0.1 & 0.07960912847730842 & -1.5598945460326698e-08\\ 50 & 0.02 & 0.0796091435489619 \\ 100 & 0.01 & 0.07960914357119007 & -1.4818785087911125e-12\\ 500 & 0.002 & 0.07960914357266899 \\ 1000 & 0.001 & 0.0796091435726711 & -1.4062824978585316e-16\\ \end{array} which confirms the overall 4th order of convergence that is inherited from the 4th order of the integration method.

Lutz Lehmann
  • 126,666