Thread: Numerical Iterations -- Predictions to achieve convergence

1. Numerical Iterations -- Predictions to achieve convergence

Using the equation:

$S = \int\limits_{a}\limits^{b}\sqrt{1 + f'(x)^2} dx$ with bounds a = 0 and b = 1; I was able to determine the arc-length of my function $f(x)$. This was done using numerical integration with enough sub-intervals to have consistent digits in the hundredths column. Knowing S I want to divide up my function $f(x)$ into equal sub-intervals of length L then I can say for the first section that:

$L = \int\limits_{0}\limits^{b}\sqrt{1 + f'(x)^2} dx$

and iterate through b to get a consistent digit in the hundredths column. This is extremely time consuming. It would be beneficial if I could integrate this equation by hand but because of the high polynomials that are incurred in y = f(x) numerical integration must be done. Is there an optimum way to guess b (similar to shooting method) to possibly get a solution in a few iterations?

Thanks

MT

2. Here's how I would approach that problem. You have a positive integrand. Therefore, $L$ is a monotone increasing function of $b$. Find a $b_{0}$ such that the integral $\int_{0}^{b_{0}}$ is larger than the desired $L$. Then utilize a bisection algorithm:

1. Split the interval $[0,b_{0}]$ in two equal pieces.
2. Evaluate the integral $\int_{0}^{b_{0}/2}.$
3. If the integral is larger than the desired L, then split the first subinterval in two. Otherwise, split the second subinterval in two.
4. Repeat steps 2 and 3 with the new interval midpoint as the upper limit of integration, continually subdividing and checking your integrals until the results differ by less than a pre-determined tolerance. Then you're done.

Example: $f(x)=x^{2}+x^{5}-3x^{7}.$ Desired $L=1.$ Tolerance in the value of L = 0.001.
Note that $\int_{0}^{2}\ge 349>2.$

1. Subdivide $[0,2]$ into $[0,1]$ and $[1,2].$
2. Evaluate $\int_{0}^{1}$. We get $2.32439.$
3. This value is greater than the desired L. Split $[0,1]$ into $[0,0.5]$ and $[0.5,1].$
2'. Evaluate $\int_{0}^{0.5}$. We get $0.578408.$
3'. Split $[0.5,1]$ into $[0.5,0.75]$ and $[0.75,1]$.
2''. Evaluate $\int_{0}^{0.75}.$ We get $0.886283.$
3''. Split $[0.75,1]$ into $[0.75,0.875]$ and $[0.875,1]$.
2'''. Evaluate $\int_{0}^{0.875}.$ We get $1.21617.$
3'''. Split $[0.75,0.875]$ into $[0.75,0.8125]$ and $[0.8125,0.875].$
2''''. Evaluate $\int_{0}^{0.8125].$ We get $0.994254.$
3''''. Split $[0.8125,0.875]$ into $[0.8125,0.84375]$ and $[0.84375,0.875]$.
2'''''. Evaluate $\int_{0}^{0.84375}.$ We get $1.08652.$
3'''''. Split $[0.8125,0.84375]$ into $[0.8125,0.828125]$ and $[0.828125,0.84375].$
2''''''. Evaluate $\int_{0}^{0.828125}.$ We get $1.03629.$
3''''''. Split $[0.8125,0.828125]$ into $[0.8125,0.820313]$ and $[0.820313,0.828125].$
2'''''''. Evaluate $\int_{0}^{0.820313}.$ We get $1.01431.$
3'''''''. Split $[0.8125,0.820313]$ into $[0.8125,0.816407]$ and $[0.816407,0.820313].$
2''''''''. Evaluate $\int_{0}^{0.816407}.$ We get $1.00406.$
3''''''''. Split $[0.8125,0.816407]$ into $[0.8125,0.814453]$ and $[0.814453,0.816407].$
2'''''''''. Evaluate $\int_{0}^{0.814453}.$ We get $0.999097.$ The error, which I've been checking as I go, is now less than the specified tolerance. We therefore pick $b=0.814453.$

You can see that our algorithm converged to an acceptable answer (meaning, within tolerance) after 9 iterations. It is also possible to specify the tolerance on the value of $b$. I'm uncomfortable with that, because you then don't know if the value of the arc length corresponding to the b you found is really all that close. With the tolerance on the L values, you know that your value of b produces an arc length that's as close as you want it to be.

All of this can be easily coded up in Mathematica or MATLAB, and I would think you'd be able to run the algorithm and get its answer in seconds at most.

3. I've got an even better idea.

Step 1. Get a hold of a copy of Mathematica (I'm not sure how you could do this in any other CAS).
Step 2. For the example I gave in Post # 2, execute the Mathematica command f[x_]:=x^2+x^5-3x^7.
Step 3. Background: We can differentiate your original equation with respect to b as follows:

$\displaystyle{\frac{dL}{db}=\sqrt{1+(f'(b))^{2}}}.$

Now, we simply invert this and view $b$ as a function of $L$. Thus,

$\displaystyle{b'(L)=\frac{1}{\sqrt{1+(f'(b(L)))^{2 }}}}.$

Viewed as a first-order DE, the initial condition becomes $b(0)=0.$ So the next step is to assign as follows:

$\text{solution}=\text{NDSolve}[\{b'[L]==\frac{1}{\sqrt{1+(f'[b[L]])^{2}}},b[0]==0\},b[L],\{L,0,5\}}]$

The result is an interpolating function.

Step 4. You can evaluate it directly like this:

$b[L]/.\text{solution}/.L\to 1$

The result is

$\{0.814814\}.$

Step 5. Checking with

$\int_{0}^{0.814814}\sqrt{1+f'[x]^{2}}\,dx//N$

yields

$1.$

I think that's pretty decent and easy convergence, don't you?

You could easily use this method to chop up a whole interval of arc into equal-sized pieces. Mathematica is very powerful in this way of symbolic computation.

5. It'll work if you can evaluate the solution b(L) at any L you want. I was just thinking about this problem. I think you could solve this in MATLAB. And if you can do it in Mathematica and MATLAB, I would think you could do it in Maple. The trick, as I said, is to be able to evaluate the numerical solution b(L) of the DE at the value of L you want, and not as an index to an array. I should also point out that for a particular problem, you might have to explicitly write out the $(f'(b(L)))^{2}.$ You follow?

Mathematica handles function definitions really well; MATLAB has an unfortunate dichotomy between its rather weak symbolic toolbox and its incredibly strong numerical, vector side. I don't know where Maple fits in that spectrum. I also don't know Maple at all, nor do I have a copy of it.

6. I have been using matlab this whole time to do the numerical integration so I guess I will stick with it for this final part. I guess the whole thing that is getting me at this point is: how do I explicitly write $(f'(b(L)))^{2}.$? I have been trying to do this using you example and I don't feel like I have a handle on the whole inversion of the function into a workable DE. I have used matlab before to solve DEs using Runge-Kutta schemes using symbolic equations. I just have to get over this hurdle i guess.

Thanks

MT

7. For my example, you'd write $(f'(b(L)))^{2}$ as follows:

First, you compute $(f'(x))^{2}.$ Well,

$f'(x)=2x+5x^{4}-21x^{6},$ and hence

$(f'(x))^{2}=(2x+5x^{4}-21x^{6})^{2}=(2x+5x^{4}-21x^{6})(2x+5x^{4}-21x^{6})$
$=4x^{2}+10x^{5}-42x^{7}+10x^{5}+25x^{8}-105x^{10}-42x^{7}-105x^{10}+441x^{12}$
$=4x^{2}+20x^{5}-84x^{7}+25x^{8}-210x^{10}+441x^{12}.$

Hence,

$(f'(b(L)))^{2}=4(b(L))^{2}+20(b(L))^{5}-84(b(L))^{7}+25(b(L))^{8}-210(b(L))^{10}+441(b(L))^{12}.$

This is why I prefer Mathematica: no need to multiply all that out. But I think here you can see how to do it.

8. I have been trying my hardest to get this to work in matlab but don't think my input equation is correct. I created a function to call into the Runge-Kutta equation solver the equation is as follows:

function dydt = f(bL)
db_dL = sqrt(1 + 4*bL^2 + 20*bL^5 - 48*bL^7 + 25*bL^8 - 210*bL^10 + 441*bL^12)

I the reason why i used bL is because I don't know a way to input $b(L)$
into matlab.

I do want to thank you for your time... I will continue to try to get matlab
to accept my input!... wish I could do it in mathematica!

MT

9. I can see at least three different problems with the code you exhibited. Perhaps fixing those problems might help MATLAB to work for you:

1. You're calling your function dydt, but you don't assign dydt in the function. You assign db_dL.
2. There's a -48*bL^7 inside the square root. It should be a -84*bL^7.
3. The derivative db/dL is actually equal to the reciprocal of the square root.

Here's some code that might work for you:

Code:
function db = arclen(L,b)
db = 1/sqrt(1 + 4*b^2 + 20*b^5 - 84*b^7 + 25*b^8 - 210*b^10 + 441*b^12)

options = odeset('RelTol',1e-4,'AbsTol',[1e-4]);
[L,b] = ode45(@arclen,[0 5],[0 0],options);
I haven't actually tested this code, since I don't have MATLAB. I'm just going by analogy with the online help.