Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
1999
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 1999

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Solving difficult integral

  • To: mathgroup at smc.vnet.net
  • Subject: [mg19094] Re: Solving difficult integral
  • From: Paul Abbott <paul at physics.uwa.edu.au>
  • Date: Thu, 5 Aug 1999 01:35:12 -0400
  • Organization: University of Western Australia
  • References: <7mek2g$9il@smc.vnet.net> <7mjr6q$f3k@smc.vnet.net> <7ni5d9$5jn@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Christian Gudrian wrote:

> Well, what am I trying to do? I want to write a C-program which does the
> following: I've got four points in R3 which define one section of a
> Catmull-Rom-Spline. With the coordinates of these four points I get a
> formula P(t) (0 <= t <= 1) which defines the curve between point two and
> thre (the points one and four are needed to calculate the tangents through
> the curve at points two and three).
>
> If I now wished (well, I certainly do!) to know the length of the path my
> point P(t) left behind at a certain point of time t I need to integrate the
> value of the velocity vector Sqrt[(dP/dt)^2] with the limits 0 and t to get
> a function s(t).
>
> But my core problem is even more complicated: I need the inverse function
> t(s) basically. What makes the whole situation even more complicated is that
> I do not know neither of the variables in advance. Right know I'm thinking
> about using MathLink to do the respective calculations during runtime when
> all things a definite. This should FullSimplify[Everthing] (:

Here is Notebook which presents one approach to your problem.

Notebook[{

Cell[CellGroupData[{
Cell["ArcLength", "Section"],

Cell[TextData[{
  "Let ",
  Cell[BoxData[
      \(TraditionalForm\`C\)]],
  " be an arc given by ",
  Cell[BoxData[
      \(TraditionalForm\`y = f(x)\)]],
  ", ",
  Cell[BoxData[
      \(TraditionalForm\`f\)]],
  " continuously differentiable in ",
  Cell[BoxData[
      \(TraditionalForm\`\([\[Alpha], \[Beta]]\)\)]],
  ". Then the length ",
  Cell[BoxData[
      \(TraditionalForm\`L\)]],
  " of ",
  Cell[BoxData[
      \(TraditionalForm\`C\)]],
  " is"
}], "Text"],

Cell[BoxData[
    \(TraditionalForm\`L = \[Integral]\_\[Alpha]\%\[Beta]\(\ at \(
                1 + \((\(f\^\[Prime]\)(x))\)\^2\)\) \
\(\(\[DifferentialD]x\)\(.\)\)\)], "DisplayFormula"],

Cell[TextData[{
  "In your case the arc length is ",
  Cell[BoxData[
      \(TraditionalForm\`s(
          t) = \[Integral]\_0\%t\(\ at \( e\ x\^4 + d\ x\^3 + c\ x\^2 + b\ x +
                  a\)\) \[DifferentialD]x\)]],
  " which can be easily computed using ",
  Cell[BoxData[
      FormBox[
        StyleBox["NIntegrate",
          "Input"], TraditionalForm]]],
  "."
}], "Text"],

Cell[BoxData[
    \(TraditionalForm\`s[a_, b_, c_, d_, e_] :=
      Function[t,
        NIntegrate[\ at \(e\ x\^4 + d\ x\^3 + c\ x\^2 + b\ x + a\), {x, 0,
            t}]]\)], "Input"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm\`\(s[1, 2, 4, 2, 1]\)[3]\)], "Input"],

Cell[BoxData[
    \(TraditionalForm\`17.114292971881074`\)], "Output"]
}, Open  ]],

Cell[TextData[{
  "However, a second approach which may be much faster is to turn this into a \

differential equation and numerically solve it over some domain (here ",
  Cell[BoxData[
      \(TraditionalForm\`t \[Element] \([0, 4]\)\)]],
  ")."
}], "Text"],

Cell[BoxData[
    \(TraditionalForm\`s[a_, b_, c_, d_, e_] := \(s[a, b, c, d, e] =
        Module[{y},
          y /. First[
              NDSolve[{\(y'\)[
                      t] == \ at \(e\ t\^4 + d\ t\^3 + c\ t\^2 + b\ t + a\),
                  y[0] == 0}, y, {t, 0, 4}]]]\)\)], "Input"],

Cell[BoxData[
    \(TraditionalForm\`\(Plot[\(s[1, 2, 4, 2, 1]\)[t], {t, 0,
          4}];\)\)], "Input"],

Cell[TextData[{
  "An easy way to construct the inverse function ",
  Cell[BoxData[
      \(TraditionalForm\`t(s)\)]],
  " is to sample the values of the plot, ",
  StyleBox["e.g.",
    FontSlant->"Italic"],
  ","
}], "Text"],

Cell[BoxData[
    \(TraditionalForm\`Cases(
      First(Plot(\(s[1, 2, 4, 2, 1]\)[t], {t, 0, 4},
          DisplayFunction -> Identity)),
      Line({x__}) -> x, \[Infinity])\)], "Input"],

Cell[TextData[{
  "which gives a list of ",
  Cell[BoxData[
      \(TraditionalForm\`{t, s}\)]],
  " values, and then use interpolation on these values after mapping ",
  Cell[BoxData[
      FormBox[
        StyleBox["Reverse",
          "Input"], TraditionalForm]]],
  " over each pair to get a list of ",
  Cell[BoxData[
      \(TraditionalForm\`{s, t}\)]],
  " values. We turn this into a procedure:"
}], "Text"],

Cell[BoxData[
    \(TraditionalForm\`t[a_, b_, c_, d_, e_] := \(t[a, b, c, d, e] =
        Interpolation(
          Reverse /@ \(Cases(
              First(Plot(Evaluate[\(s[a, b, c, d, e]\)[t]], {t, 0, 4},
                  DisplayFunction -> Identity)),
              Line({x__}) -> x, \[Infinity])\))\)\)], "Input",
  AspectRatioFixed->True,
  CellTags->"InverseFunction"],

Cell[TextData[{
  "For the above curve, here is the ",
  Cell[BoxData[
      \(TraditionalForm\`t\)]],
  "-value corresponding to ",
  Cell[BoxData[
      \(TraditionalForm\`s = 21\)]],
  "."
}], "Text"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm\`\(t[1, 2, 4, 2, 1]\)[21]\)], "Input",
  AspectRatioFixed->True],

Cell[BoxData[
    \(TraditionalForm\`3.271339339542367`\)], "Output"]
}, Open  ]],

Cell[TextData[{
  "and here is a plot of ",
  Cell[BoxData[
      \(TraditionalForm\`t(s)\)]],
  "."
}], "Text"],

Cell[BoxData[
    \(TraditionalForm\`\(Plot[\(t[1, 2, 4, 2, 1]\)[s], {s, 0,
          34}];\)\)], "Input"]
}, Open  ]]
}
]


Cheers,
    Paul
____________________________________________________________________
Paul Abbott                                   Phone: +61-8-9380-2734
Department of Physics                           Fax: +61-8-9380-1014
The University of Western Australia
Nedlands WA  6907                     mailto:paul at physics.uwa.edu.au
AUSTRALIA                            http://physics.uwa.edu.au/~paul

            God IS a weakly left-handed dice player
____________________________________________________________________




  • Prev by Date: Re: equaltity of lists
  • Next by Date: Re: Biased Random[Integer]?
  • Previous by thread: Re: compiling problem
  • Next by thread: Re: How to prevent a window from being resized?