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 ____________________________________________________________________