From 32376186d729de62094dcf0131ebd00854e08618 Mon Sep 17 00:00:00 2001 From: Thomas Baruchel Date: Tue, 4 Oct 2022 18:30:41 +0200 Subject: [PATCH] Update --- thiele.mac | 64 ------------------------------------------------------ 1 file changed, 64 deletions(-) delete mode 100644 thiele.mac diff --git a/thiele.mac b/thiele.mac deleted file mode 100644 index 6124647..0000000 --- a/thiele.mac +++ /dev/null @@ -1,64 +0,0 @@ -/* -Copyright (c) 2020 Thomas Baruchel - -Permission is hereby granted, free of charge, to any person obtaining a copy of -this software and associated documentation files (the "Software"), to deal in -the Software without restriction, including without limitation the rights to -use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies -of the Software, and to permit persons to whom the Software is furnished to do -so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - -/* - - In the following example, we use Thiele's interpolation formula - for computing the inverse cos function cos^(-1) from 32 points - given by the cos function itself. - - We create two vectors: - * the vector 'a' contains regularly spaced points in [0..1.55] - * the vector 'b' contains cos(x) for each x in vector 'a' - - We set the keepfloat flag to true in order to prevent Maxima - to replace the values with rational approximations. - - Then we compute the inverse function by using thiele(b,a); - we evaluate this interpolation at 0.5 expecting some value - close to pi/3. - -(%i18) a:makelist(0.05*i,i,0,31)$ -(%i19) b:makelist(float(cos(0.05*i)),i,0,31)$ -(%i20) keepfloat:true$ -(%i21) subst(x=0.5, thiele(b,a))*3; -(%o21) 3.14159265357928 - -*/ - - -/* Thiele's interpolation formula */ -thiele(u, v) := block([rho:makelist( - makelist(v[i], length(v)-i+1), - i, length(v)), a:0], - for i:1 thru length(rho)-1 - do rho[i][2]: (u[i]-u[i+1]) / (rho[i][1] - rho[i+1][1]), - for i:3 thru length(rho) - do (for j:1 thru length(rho)-i+1 - do rho[j][i]: (u[j]-u[j+i-1]) - / (rho[j][i-1]-rho[j+1][i-1]) - + rho[j+1][i-2]), - rho: rho[1], - for i:length(rho) thru 3 step -1 - do a: ratsimp(( 'x - u[i-1])/(rho[i]-rho[i-2]+a)), - ratsimp( v[1] + ( 'x - u[1] ) / (rho[2] + a) ))$