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[
" be an arc given by ",
Cell[BoxData[
", ",
Cell[BoxData[
" continuously differentiable in ",
Cell[BoxData[
". Then the length ",
Cell[BoxData[
" of ",
Cell[BoxData[
" is"
}], "Text"],

Cell[BoxData[
1 + \((\(f\^\[Prime]\)(x))\)\^2\)\) \
\(\(\[DifferentialD]x\)\(.\)\)\)], "DisplayFormula"],

Cell[TextData[{
"In your case the arc length is ",
Cell[BoxData[
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",
"."
}], "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[
}, 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[
")."
}], "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[
" is to sample the values of the plot, ",
StyleBox["e.g.",
FontSlant->"Italic"],
","
}], "Text"],

Cell[BoxData[
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[
" values, and then use interpolation on these values after mapping ",
Cell[BoxData[
FormBox[
StyleBox["Reverse",
" over each pair to get a list of ",
Cell[BoxData[
" 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[
"-value corresponding to ",
Cell[BoxData[
"."
}], "Text"],

Cell[CellGroupData[{

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

Cell[BoxData[
}, Open  ]],

Cell[TextData[{
"and here is a plot of ",
Cell[BoxData[
"."
}], "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?