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
____________________________________________________________________