Results 1 to 9 of 9

Math Help - Numerical Iterations -- Predictions to achieve convergence

  1. #1
    Newbie
    Joined
    Nov 2009
    Posts
    19

    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
    Follow Math Help Forum on Facebook and Google+

  2. #2
    A Plied Mathematician
    Joined
    Jun 2010
    From
    CT, USA
    Posts
    6,318
    Thanks
    4
    Awards
    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.
    Follow Math Help Forum on Facebook and Google+

  3. #3
    A Plied Mathematician
    Joined
    Jun 2010
    From
    CT, USA
    Posts
    6,318
    Thanks
    4
    Awards
    2
    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.
    Follow Math Help Forum on Facebook and Google+

  4. #4
    Newbie
    Joined
    Nov 2009
    Posts
    19
    I don't have access to mathematica but I do have access to maple 14 would that work?
    Follow Math Help Forum on Facebook and Google+

  5. #5
    A Plied Mathematician
    Joined
    Jun 2010
    From
    CT, USA
    Posts
    6,318
    Thanks
    4
    Awards
    2
    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.
    Follow Math Help Forum on Facebook and Google+

  6. #6
    Newbie
    Joined
    Nov 2009
    Posts
    19
    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
    Follow Math Help Forum on Facebook and Google+

  7. #7
    A Plied Mathematician
    Joined
    Jun 2010
    From
    CT, USA
    Posts
    6,318
    Thanks
    4
    Awards
    2
    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.
    Follow Math Help Forum on Facebook and Google+

  8. #8
    Newbie
    Joined
    Nov 2009
    Posts
    19
    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
    Follow Math Help Forum on Facebook and Google+

  9. #9
    A Plied Mathematician
    Joined
    Jun 2010
    From
    CT, USA
    Posts
    6,318
    Thanks
    4
    Awards
    2
    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.
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Number of Iterations
    Posted in the Discrete Math Forum
    Replies: 6
    Last Post: September 16th 2010, 06:54 AM
  2. Iterations ??
    Posted in the Calculus Forum
    Replies: 1
    Last Post: June 10th 2010, 02:34 AM
  3. Replies: 1
    Last Post: May 14th 2010, 09:04 AM
  4. Achieve Linearity
    Posted in the Advanced Statistics Forum
    Replies: 1
    Last Post: December 1st 2009, 11:24 PM
  5. Replies: 4
    Last Post: September 18th 2009, 02:16 PM

Search Tags


/mathhelpforum @mathhelpforum