27.4.3.1 Choosing a Boundary Value

If we return to our problem, we have a boundary condition at r = 0 and one boundary condition at r = R. We now decide from which boundary we want to use our numerical schemes. If you look at Eq. 27.51, you can see that the independent variable r occurs in the denominator. This is a very good indication that choosing the boundary for r = 0 is not a good idea, as this would immediately produce a division by zero. Obviously, if this is the first value we calculate in our numerical scheme, it will render all subsequent calculations impossible. However, if this value occurs as the last value in our calculation, it may be significantly less problematic, as no subsequent calculations would have to deal with an infinite value. We have already encountered such a scenario when solving the contour of sessile and pendant drops numerically. Refer to section 24.4.2.3 for details.

These cases occur pretty often, especially when working in spherical or cylindrical coordinates. It is therefore, a good idea to check precisely which boundary is the best to start the calculation from. Obviously, choosing the boundary will not affect the result in terms of precision. Proceeding end-to-front should produce the same results as proceeding front-to-end. However, in many cases, one of these directions is not possible for reasons such as the one we just encountered.

Therefore we opt for r = R as the boundary condition from which we want to proceed. However, we still have not overcome the problem of the mixed boundary condition.

27.4.3.2 Rationale

In our case, we require a boundary value for u (r = R). The question is, how can we derive this value from the equations we have? The answer is simple: we cannot. We can only guess the value. Using a guessed value, we will make a “shot” at the problem and see, which value of u (r = 0) our guess for u (r = R) ends up producing. If this value is far from the correct value, we take a second guess checking if this guess produces a better result.

This is the reason why this method is referred to as shooting. It is pretty much an analogy to shooting arrows. If your first arrow is off to the left of the target, you would aim your second shot a bit further to the right. If the second shot is too far right of the target, you would aim your next shot a bit further to the left again. We can do this until we end up with a guess for u (r = R), which produces a value for u (r = 0) which is reasonably close to the value we expect. In our case, this value would be 0.

27.4.3.3 Implementation

We now use the method of shooting for solving Eq. 27.50 and Eq. 27.51. The implementation is shown in listing 27.9. In line 1 we restart the Maple server and include Core.txt and CoreFunctions.txt the latter of which contains the function RungeKutta4ODEs. The next line contains a directive that tells Maple to display an array with all elements. Usually, Maple will only display an array in abbreviated form once it becomes too large. However, it is easier to see the numerical results if we can see the whole matrix.

1 restart:read "Core.txt":read "CoreFunctions.txt":
2 interface(rtablesize=infinity):
3
4 guess1:=RungeKutta4ODEs([(r,v,u)->u,(r,v,u)->-10E3-u/r],[v,u],0.025,0,[0,1],11);
5 guess2:=RungeKutta4ODEs([(r,v,u)->u,(r,v,u)->-10E3-u/r],[v,u],0.025,0,[0,100],11);
6 guess3:=RungeKutta4ODEs([(r,v,u)->u,(r,v,u)->-10E3-u/r],[v,u],0.025,0,[0,−100],11);
7 guess4:=RungeKutta4ODEs([(r,v,u)->u,(r,v,u)->-10E3-u/r],[v,u],0.025,0,[0,−125],11);
8
9 display(quickplot([r->-1.5625*((r/0.025)^2-1)],0..0.025,typeset("r",[mm]),typeset(v[x]^empty,[mm/s])),pointplot(extracttoplot(guess3,[1,2]),symbolsize=7,color=black,symbol=solidcircle));
10 display(quickplot([r->-1.5625*((r/0.025)^2-1)],0..0.025,typeset("r",[mm]),typeset(v[x]^empty,[mm/s])),pointplot(extracttoplot(guess4,[1,2]),symbolsize=7,color=black,symbol=solidcircle));

Line 4 shows our first guess. We call the function RungeKutta4ODEs with Eq. 27.50 and Eq. 27.51 (in the correct order) and the two dependent variables v and u. We indicate the range to go from r = 0.025mm to r = 0mm, i.e., back-to-front in order to avoid the problem of the division by zero. Just to illustrate, you may want to change the order, in which case Maple will abort the method immediately with a division by zero error. The next argument is the critical point: the initial values. We know the exact initial value for v (R) = 0 (which is zero due to the no-slip boundary condition), but we do not know the boundary condition for u (R). Here, the first guess is 1 which is as good as any guess. The last parameter indicates that we want a total of 11 points to be calculated, which is not a lot, but sufficient.

If you run line 4, you will see that the first guess is horribly wrong. The array displayed shows us that the value the algorithm finds for the second boundary is u (0) = 8. This value results from the division by zero and would not be problematic. However, if you look at the second-to-last value you will see that it is u (0.0025) = 1247.5.

This value is definitely not approaching 0 and, in fact if we look at the calculated values of u (r), we can see that they are growing exponentially. This guess is therefore definitely wrong.

Our second guess (see line 5) makes things even worse as u (0.0025) = 2237.5 and the function is growing even more rapidly. Therefore it seems reasonable to assume a negative value for u (R). In fact, as we know the analytical solution given by Eq. 27.49, we know that this value must be negative. We even know its exact value to be

$\begin{array}{ll}u\left(r\right)\hfill & =\frac{dv}{dr}=-1.5625\mathrm{mm}{\mathrm{s}}^{-1}·2\frac{r}{R^2}\hfill \\ u\left(R\right)\hfill & =-3.125\mathrm{mm}{\mathrm{s}}^{-1}\frac{1}{25\mathrm{\mu m}}=-125{\mathrm{s}}^{-1}\hfill \end{array}$

However, in most applications, we do obviously not have access to the correct solutions so shooting is the only option. Our third guess (see line 6) is definitely better. However the function u (r) still grows where it should actually converge to zero. If we plot the output of this guess together with the analytical solution, we obtain Fig. 27.14a. As you can see, the solution is off from the correct result. If we run the last guess (see line 7), we find that the function u (r) now converges to zero. If we plot the results from this guess, we obtain Fig. 27.14b. As you can see, this is the correct result. The plots are produced from line 9 and line 10.

27.4.3.4 Pendant Drop Revisited

Introduction. We will introduce a second example of a shooting method using a system of ODEs that we have derived during the study of the pendant drop (see section 24.5). We have used a numerical scheme (see listing 24.2) to solve the rather complicated ODE given by Eq. 24.10, which describes the shape of the drop. For the approach of del Río and Neumann, we required Eq. 24.22, Eq. 24.23, Eq. 24.24, and Eq. 24.16 given by

$\begin{array}{lll}\frac{d\theta }{ds}\hfill & =\hfill & -\frac{sin\theta }{x}+\frac{\mathit{∆\rho gz}}{\mathrm{\gamma }}+\frac{2}{r_{\text{max}}}\hfill \\ \frac{dx}{ds}\hfill & =\hfill & cos\theta \hfill \\ \frac{dz}{ds}\hfill & =\hfill & sin\theta \hfill \\ \frac{dV}{ds}\hfill & =\hfill & \pi x^{2}sin\theta \hfill \\ \frac{drmax}{ds}\hfill & =\hfill & 0\hfill \\ \hfill & \hfill & \hfill \end{array}$

At the time, we noted that using the path variable s is actually impractical because this value tells us little about our problem, and we do not require it as a result. Also it would be desirable to rewrite these equations as a function of the volume V of the drop because this is the actual independent variable. We can rewrite these equations using

$ds=\frac{dV}{p{x}^{2}sin\theta }$

in which case we find the following system of ODEs

$\frac{d\theta }{dV}=-\frac{1}{x^3\mathrm{\pi }}+\frac{\mathrm{\Delta }\rho g}{\mathrm{\gamma }}\frac{z}{{\mathrm{\pi }\mathrm{x}}^{2}sin\theta }+\frac{2}{r_{\text{max}}\mathrm{\pi }{x}^{2}sin\theta }$

$\frac{dx}{dV}=\frac{cos\theta }{\mathrm{\pi }{x}^{2}sin\theta }=\frac{1}{\mathrm{\pi }{x}^{2}tan\theta }$

$\frac{dz}{ds}=\frac{sin\theta }{\mathrm{\pi }{x}^{2}sin\theta }=\frac{1}{\mathrm{\pi }{x}^{2}sin\theta }$

$\frac{d{r}_{\text{max}}}{dV}=0$

where we have dropped the ODE for the path variable s, as we are not interested in this value.

Boundary Values. Obviously, our problem has a number of boundary values given by

$\begin{array}{llll}\theta \left(0\right)\hfill & =0,\hfill & \lambda \left(V_{\mathit{max}}\right)\hfill & =\to \hfill \\ x\left(0\right)\hfill & =0,\hfill & x\left(V_{\mathit{max}}\right)\hfill & =r_0\hfill \\ z\left(0\right)\hfill & =0,\hfill & z\left(V_{\mathit{max}}\right)\hfill & =\to \hfill \\ r_{\text{max}}\left(0\right)\hfill & =0,\hfill & r_{\text{max}}\left(V_{\mathit{max}}\right)\hfill & =\to \hfill \end{array}$

where V_max (the volume of the drop) and r₀ (the diameter of the capillary to which it is suspended) are given and known values. Obviously, it would be ideal to start from the lower boundary, i.e., going from V = 0 to V = V_max, as we know almost all boundary conditions and merely have to determine r_max. However, if we look at our ODEs, we can see that we find the value of x in several denominators therefore we would be faced with division by zero errors.

Therefore we will have to proceed from back-to-front starting from V = V_max. Unfortunately, there are quite a number of values which we do not know and we have to shoot for three unknown boundary conditions at V = V_max.

Implementation. Listing 27.10 shows the implementation of this shooting approach. It starts by implementing Eq. 27.54, Eq. 27.55, Eq. 27.56, and Eq. 27.57 (see line 1). The term $\frac{\mathrm{\Delta }\rho g}{\mathrm{\gamma }}$ is converted to the value - 0.13625 using the values for water, keeping in mind that the gravity for pendant drop is negative. We then apply the function RungeKutta40DEs from CoreFunctions.txt in line 6. We pass the four ODEs and the range for V. Here we use a water droplet with a total volume of 18.24 µl. This example is taken from listing 24.2 - it is the largest droplet studied in this listing. We use a total of 201 data points and pass the initial values we know which is only x (V_max) = r₀ = 0.5. This value is also taken from the example in listing 24.2.

1 dphidV:=(V,phi,x,z,rmax)->-1/(x^3*Pi)-0.13625*z/(Pi*x^2*sin(phi))+2/(rmax*Pi*x^2*sin(phi));
2 dxdV:=(V,phi,x,z,rmax)->1/(Pi*x^2*tan(phi));
3 dzdV:=(V,phi,x,z,rmax)->1/(Pi*x^2);
4 drmaxdV:=(V,phi,x,z,rmax)->0;
5
6 sol:=RungeKutta40DEs([dphidV,dxdV,dzdV,drmaxdV],[phi,x,z,rmax ],18.24,0,[1.445,0.5,4.315,1.4356],201):sol [201];
7
8 quickplot([extracttoplot(sol,[3,4])],0..5,"x [mm]","z [mm]");

You will note that the output is written in a way that only the last row of the array created will be displayed (see line 6). The reason for this is that we are only interested in the values of this last row, which gives the values of the boundary for V = 0. These values must be close to the known values for the boundary. If the shooting was correct, we should see a similar output here. For the boundary values shown in listing 27.10, line 6 we obtain the following values at the lower boundary

$\begin{array}{l}\theta \left(0\right)=0.022˜0\\ x\left(0\right)=0.036˜0\\ z\left(0\right)=-0.018˜0\end{array}$

which are all sufficiently close to the desired values. Fig. 27.15 compares the values produced by our shooting method with the values we obtained previously in listing 24.2. As you can see, the results are very close.

Speeding Up the Shooting Process. It may take a significant amount of time to apply the shooting method to a complicated case such as the one shown here. Shooting three unknown boundary values obviously is a purely trial-and-error search for a good set of values. However, this process can be automated because an algorithm can try a wide range of values and determine the ones that yield the best result. This, in fact, is a strength of shooting. The quality of the result can be assessed by simply checking how close the produced values are to the desired ones. Obviously, in a case such as the one described here, searching manually is a very tedious process.

27.4.3.5 Summary

As we have seen in this section, the method of shooting is one of the easiest methods to implement whenever an ODE or system of ODEs with mixed boundary conditions is encountered. However, it may require a significant amount of trials until a good solution is obtained. In many application scenarios, it is reasonable to assume that a vague idea about the size or at least the sign of the unknown second boundary condition is available. This speeds up the solution process significantly. We have seen one example derived from the Navier-Stokes equation in cylindrical coordinates that we were able to solve sufficiently exact using a shooting approach to derive the unknown boundary condition.

27.4.4 Function Approximation Methods

27.4.4.1 Introduction

The method of shooting is not the only method we can use when we have mixed boundary conditions. One very important set of methods is referred to as function approximation methods. As the name suggests, we use a function that we assume can approximate the solution to the ODE. We usually refer to this function as the trial function. There are many different types of trial functions one can use, but they will contain some unknown parameters P₀, P_1, . . . , P_n which need to be determined. The trial function is usually denoted T (x, P₀, P1, … , P_n) where x is the independent variable and P₀, P₁, . . . , P_n are the parameters.

In the next step, we would substitute the trial function into the ODE. If we are very lucky, we will find that the trial function solves the ODE. This would make the trial function a solution immediately. As can be imagined, this case rarely happens in practical applications. We would rather be left with an equation that has n + 1 unknown and, which (being a single equation) will therefore not yield a solution at all.

Collocation Methods. This is where we will use a collocation method. Instead of demanding that the substituted trial function will solve the original ODE on the whole domain of x, we simply demand that it satisfies it at n + 1 discreet values x_i. Evaluating the substituted trial function at these values will give a set of n + 1 equations, which we have to solve for n + 1 unknowns, i.e., the parameters P_i. At each discreet value x_i the ODE would then be satisfied. This is why these approaches are referred to as collocation methods. At n + 1 distinct collocation points, our solution will be exact. Obviously, the higher the number of parameters used, the more collocation points we will have and the more exact the solution.

Example. In order to illustrate this example, we will use the Navier-Stokes equation in cylindrical coordinates Eq. 27.45 again, which is given by

$\frac{1}{r}\frac{dv}{dr}+\frac{d^{2}v}{d{r}^2}=\frac{1}{\eta }\frac{dp}{dz}$

Obviously, we know the analytical solution to this ODE, which is given by Eq. 16.20. As we will see, we can also derive this solution using a collocation approach.

27.4.4.2 Trial Function

A very common type of trial functions are polynomial series for one simple reason: they are linear and can be linearly differentiated. As an example, consider the following trial function

$T\left(x,p_0,p_1,\dots ,p_n\right)={\displaystyle \underset{i=0}{\overset{n}{\Sigma }}{p}_i{f}_i\left(x\right)}$

where f_i (r) are functions which are chosen because we assume that the solution may be well approximated by them. Given our ODE, we may assume that a suitable series of functions would be

$T\left(r,p_0,p_1,\dots ,p_n\right)={\displaystyle \underset{i=0}{\overset{n}{\Sigma }}{p}_i}{r}^i$

If we differentiate this trial function with respect to the independent variable, we obtain the general form

$\frac{dT}{dx}={\displaystyle \underset{i=0}{\overset{n}{\Sigma }}{p}_i\frac{d{f}_i}{dx}\left(x\right)}$

This means that irrespective of the type of function f_i, we obtain again a polynomial but this time of the

derivatives $\frac{d{f}_i}{dx}$. For our assumed solution we would obtain

$\frac{dT}{dx}\left(r,p_0,p_1,\dots ,p_n\right)={\displaystyle \underset{i=1}{\overset{n}{\Sigma }}{p}_{i}i{r}^{i-1}}$

27.4.4.3 Residual Function

Obviously, the trial function must be applied to the ODE in order to test whether or not it is a good approximation for a solution. The resulting function is referred to as residual function ? (x, p_0,p₁, . . . , p_n). If we apply our trial function Eq. 27.60 to our ODE Eq. 27.58 we obtain

$\frac{1}{r}{\displaystyle \underset{i=1}{\overset{n}{\Sigma }}{p}_{i}i{r}^{i-1}+}{\displaystyle \underset{i=2}{\overset{n}{\Sigma }}{p}_{i}i\left(i-1\right)r^{i-2}}=\frac{1}{\eta }\frac{dp}{dz}$

As you can see, we are now left with an ordinary equation as the differentials have disappeared. Unfortunately, we have the unknowns p_i for which we do not have values at the moment. However, we are now free to evaluate this function for arbitrary values r_i. If we do this for n values of r_i, this leaves us with n residual equations. We now demand that at each of these points, Eq. 27.61 is fulfilled, i.e., the left-hand and the right-hand sides match. This gives us the required number of equations to solve for all unknowns P_i.

27.4.4.4 Adding Boundary Conditions

One point to keep in mind is that we do not require n equations because we have boundary conditions. The no-slip boundary condition allows us to set our trial function to zero for r = R

$T\left(R,p_0,p_1,\dots ,p_n\right)={\displaystyle \underset{i=0}{\overset{n}{\Sigma }}{p}_i{R}^i=0}$

The symmetry condition allows us to set the derivative $\frac{dT}{dr}=0$ to zero at r = 0

$\frac{dT}{dx}\left(0,p_0,p_1,\dots p_n\right)={\displaystyle \underset{i=1}{\overset{n}{\Sigma }}{p}_{i}i}$

This means that if we have a total of n + 1,we only require n - 1 additional equations as we already have two equations derived from the boundary conditions.

27.4.4.5 Example

In the following, we will illustrate the method using a trial function with n = 3 given by

$\begin{array}{ll}T\left(r,p_0,p_1,p_2,p_3\right)\hfill & =p_0+p_{1}r+p_2{r}^2+p_3{r}^3\hfill \\ \frac{dT}{dr}\left(r,p_0,p_1,p_2,p_3\right)\hfill & =p_1+2p_{2}r+3p_3{r}^2\hfill \\ \frac{d^{2}T}{d{r}^2}\left(r,p_0,p_1,p_2,p_3\right)\hfill & =2p_2+6p_3{r}^2\hfill \end{array}$

which we resubstitute into Eq. 27.58 to find the residual function

$\frac{1}{r}\left(p_1+2p_{2}r+3p_3{r}^2\right)+2p_2+6p_{3}r=\frac{1}{\lambda }\frac{dp}{dz}$

We now have the two boundary conditions, which give rise to the following two equations

$p_0+p_{1}R+p_2{R}^2+p_3{R}^3=0\mathrm{no}-\mathrm{slip}$

$p_1+2p_{2}0+3p_3{0}^2=0\mathrm{symmetry}$

We therefore only require two additional equations. We choose randomly, r₁= $\frac{1}{3}$R and r₂ = $\frac{2}{3}$ in which case we find from Eq. 27.63 the two equations

$\frac{3}{R}\left(p_1+2p_2+\frac{1}{3}R+3p_3{\left(\frac{1}{3}R\right)}^2\right)+2p_2+6p_3\frac{1}{3}R=\frac{1}{\eta }\frac{dp}{dz}$

$\frac{3}{2R}\left(p_1+2p_2+\frac{2}{3}R+3p_3{\left(\frac{2}{3}R\right)}^2\right)+2p_2+6p_3\frac{2}{3}R=\frac{1}{\eta }\frac{dp}{dz}$

which can be simplified to

$\frac{3}{R}{p}_1+4p_2+3R{p}_3=\frac{1}{\eta }\frac{dp}{dz}$

$\frac{3}{2R}{p}_1+4p_2+6R{p}_3=\frac{1}{\eta }\frac{dp}{dz}$

In general, the four equations Eq. 27.64, Eq. 27.65, Eq. 27.66, and Eq. 27.67 make up a linear system of equations which can be solved by any method discussed in section 25.2 or section 25.3. However, in this case, the system is so simple that it can be solved by hand. From Eq. 27.65 we can see that P₁= 0 in which case we are left with the following three equations

$\begin{array}{lllllll}{p}_0\hfill & +\hfill & R^2{p}_2\hfill & +\hfill & R^3{p}_3\hfill & =\hfill & 0\hfill \end{array}$

$\begin{array}{lllllll}\hfill & +\hfill & 4p_2\hfill & +\hfill & 3R{p}_3\hfill & =\hfill & \frac{1}{\eta }\frac{dp}{dz}\hfill \end{array}$

$\begin{array}{lllllll}\hfill & +\hfill & 4p_2\hfill & +\hfill & 6R{p}_3\hfill & =\hfill & \frac{1}{\eta }\frac{dp}{dz}\hfill \end{array}$

Subtracting Eq. 27.69 from Eq. 27.70 yields

$3R{p}_3=0\to p_3=0$

in which case we find from Eq. 27.69 that

$p_2=\frac{1}{4\eta }\frac{dp}{dz}$

and from Eq. 27.68 that

$p_0=-R^2{p}_2=-\frac{R^2}{4\eta }\frac{dp}{dz}$

At this point, all unknowns have been determined and we can note down the final solution according to Eq. 27.62

$\begin{array}{ll}v\left(r\right)\hfill & =T\left(r,p_0,p_1,p_2,p_3\right)\hfill \\ \hfill & =p_0+p_{1}r+p_2{r}^2+p_3{r}^3\hfill \\ \hfill & =\frac{1}{4\eta }\frac{dp}{dz}\left(r^2-R^2\right)\hfill \\ \hfill & =\frac{R^2}{4\eta }\frac{dp}{dz}\left({\left(\frac{r}{R}\right)}^2-1\right)\hfill \end{array}$

which obviously is the correct solution (see Eq. 16.20).

Assessment. This result is somewhat astonishing as we did (in the classical sense) not even use numerical methods to solve it. This is due to the nature of the trial function we used. We picked the correct type of collocation function which allowed us to simply determine the missing coefficients. Obviously, if we had chosen a different collocation function we would have been left with somewhat different results, which would have required a numerical derivation.

27.4.4.6 Summary

In this section we have discussed function approximation methods and how they are used to solve ODEs with boundary conditions. Approximation methods are a large class of methods that are used in many aspects of numerical modeling. We have briefly outlined their use using the Navier-Stokes equation in cylindrical coordinates as exemplary ODE for which we were able to derive the correct solution.