Thread: Numerical Iterations -- Predictions to achieve convergence

Using the equation:

with bounds a = 0 and b = 1; I was able to determine the arc-length of my function . 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 into equal sub-intervals of length L then I can say for the first section that:



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

Here's how I would approach that problem. You have a positive integrand. Therefore, is a monotone increasing function of . Find a such that the integral is larger than the desired . Then utilize a bisection algorithm:

1. Split the interval in two equal pieces.
2. Evaluate the integral
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: Desired Tolerance in the value of L = 0.001.
Note that

1. Subdivide into and
2. Evaluate . We get
3. This value is greater than the desired L. Split into and
2'. Evaluate . We get
3'. Split into and .
2''. Evaluate We get
3''. Split into and .
2'''. Evaluate We get
3'''. Split into and
2''''. Evaluate We get
3''''. Split into and .
2'''''. Evaluate We get
3'''''. Split into and
2''''''. Evaluate We get
3''''''. Split into and
2'''''''. Evaluate We get
3'''''''. Split into and
2''''''''. Evaluate We get
3''''''''. Split into and
2'''''''''. Evaluate We get The error, which I've been checking as I go, is now less than the specified tolerance. We therefore pick

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 . 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.

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:



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



Viewed as a first-order DE, the initial condition becomes So the next step is to assign as follows:



The result is an interpolating function.

Step 4. You can evaluate it directly like this:



The result is



Step 5. Checking with



yields



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.

I don't have access to mathematica but I do have access to maple 14 would that work?

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 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.

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 ? 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

For my example, you'd write as follows:

First, you compute Well,

and hence





Hence,



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

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
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

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.