Parabolic interpolation with Maple

154 Views Asked by At

Given the points of coordinates $(k, 2)$, $(\ell, 4)$, $(m, 8)$, I'm looking for:

  • the parabolic polynomial $f(x)=a x^2 + b x + c$ that goes by these points
  • and more specifically the value $n$ such that $f(n) = 0$

Of course, I know how to do this manually (simple exercise), but how to find this via Maple (or any online tool, not requiring to install a software locally)?

1

There are 1 best solutions below

2
On BEST ANSWER
restart;

f:=y=a*x^2+b*x+c:

pts:=[[k,2],[l,4],[m,8]]:

eqs:=map[2](eval,f,map[2](Equate,[x,y],pts)):

lprint(eqs);
  [2 = a*k^2+b*k+c, 4 = a*l^2+b*l+c, 8 = a*m^2+b*m+c]

sol:=simplify(eliminate(eqs,{a,b,c})[1]):

lprint(sol);
  {a = (4*k-6*l+2*m)/((l-m)*(k-m)*(k-l)),
   b = (-4*k^2+6*l^2-2*m^2)/((l-m)*(k-m)*(k-l)),
   c = ((8*l-4*m)*k^2+(-8*l^2+4*m^2)*k+2*m*l*(l-m))/((l-m)*(k-m)*(k-l))}

spec:=simplify([solve(eval(eval(f,sol),y=0),{x})]):

lprint(spec);
  [{x = (2*k^2-3*l^2+m^2+4^(1/2)*(((9/2)*l^2+(-6*k-3*m)*l+k^2+4*k*m-(1/2)*m^2)*((1/2)*l^2+(-2*k+m)*l+k^2-(1/2)*m^2))^(1/2))/(4*k-6*l+2*m)},
   {x = (2*k^2-3*l^2+m^2-4^(1/2)*(((9/2)*l^2+(-6*k-3*m)*l+k^2+4*k*m-(1/2)*m^2)*((1/2)*l^2+(-2*k+m)*l+k^2-(1/2)*m^2))^(1/2))/(4*k-6*l+2*m)}]

And now for a numeric example, using those solutions.

# Numeric example
vals:=[k=6,l=3,m=5]:

eval(pts,vals);

                [[6, 2], [3, 4], [5, 8]]

eval(sol,vals);

               /    -8      70         \ 
              { a = --, b = --, c = -42 }
               \    3       3          / 

eval(eval(y,f),eval(sol,vals));

                     8  2   70       
                   - - x  + -- x - 42
                     3      3        

xspec:=eval(spec,vals):

lprint(xspec);
  [{x = 35/8+(1/16)*4^(1/2)*217^(1/2)},
   {x = 35/8-(1/16)*4^(1/2)*217^(1/2)}]

evalf(xspec);

         [{x = 6.216364982}, {x = 2.533635018}]

And we could plot it.

plots:-display(
  plot(eval(eval(y,f),eval(sol,vals)),x=0..8),
  plot(eval(pts,vals),
       style=point,color=magenta,
       symbol=solidcircle,symbolsize=15),
  plot([[eval(x,xspec[1]),0],[eval(x,xspec[2]),0]],
       style=point,color=blue,
       symbol=solidcircle,symbolsize=15),
view=-1..10);

enter image description here

[edited below] I was asked to clarify some of the computations. The code starts with the following (and now I will not suppress output),

restart;

f:=y=a*x^2+b*x+c;

                               2
                   f := y = a x  + b x + c

pts:=[[k,2],[l,4],[m,8]];

             pts := [[k, 2], [l, 4], [m, 8]]

I want to substitute in values for x and y in that formula for f, using each of the points in the list pts above.

Let's first look at this preliminary part. For each of the three "points", such as the first [k,2], it constructs a list of the equations, such as [x=k, y=2]. Here's how just the first point could be handled, alone.

 Equate( [x,y], [k,2] );

                   [x = k, y = 2]

But I want a one-liner that produces all three at once. I can use the map or seq commands for doing that. Each of the following constructs the same thing (and my original code happened to use the first of these approaches).

map[2]( Equate,[x,y], pts );

         [[x = k, y = 2], [x = l, y = 4], [x = m, y = 8]]

map( p -> Equate([x,y],p), pts);

         [[x = k, y = 2], [x = l, y = 4], [x = m, y = 8]]

[ seq( Equate([x,y],p) , p in pts )];

         [[x = k, y = 2], [x = l, y = 4], [x = m, y = 8]]

The above list of lists of equations can be used to substitute into f for provide values for x and y.

Here's how it could be done for just the first point. (I only need to do the second of these two commands. The first is shown just to illustrate what is being passed as the 2nd argument to the eval command.

Equate( [x,y], [k,2] );

                   [x = k, y = 2]

eval( f, Equate( [x,y], [k,2] ) );

                       2
                2 = a k  + b k + c

But I wanted a one-liner that produces a list of three equations, where each is the result of subsituting for x and y into f based on each of the three points.

Here again I'll show three ways to construct the thing I want in a one-liner. I happened to use the first of these originally, just because I think it's terse.

map[2]( eval, f, map[2](Equate,[x,y],pts) );

           2                   2                   2
   [2 = a k  + b k + c, 4 = a l  + b l + c, 8 = a m  + b m + c]

map( e -> eval(f,e), map[2](Equate,[x,y],pts) );

           2                   2                   2
   [2 = a k  + b k + c, 4 = a l  + b l + c, 8 = a m  + b m + c]

[ seq( eval(f,e), e in  map[2](Equate,[x,y],pts) ) ];

           2                   2                   2
   [2 = a k  + b k + c, 4 = a l  + b l + c, 8 = a m  + b m + c]

I hope that helps explain the initial construction of the three equations, which I then solved using the eliminate command.