f# - 获取拉格朗日插值的多项式表示

Question

如何表示不完整的数学函数？

我需要做类似(x - 常量) 的事情，然后

(x - constant)*(x - another) => (x^2 - x * constant - x * another + constant * another)

等等。

我正在尝试编写一个程序来进行拉格朗日插值(为某些点找到一个函数)所以我需要从一组已知值中创建一个我可以看到的函数(打印或其他东西)。抱歉，如果令人困惑。

Answer 1

如果您想实现所讨论的拉格朗日插值here

获取插入值的函数:

那么这是直接翻译成 F#:

let LagrangeInterpol (points : (Double*Double)[]) x =
    let indizes = [0..points.Length-1]
    let p j = 
        indizes
        |> List.map (fun k -> 
            if k <> j 
            then (x - fst points.[k]) 
                  / (fst points.[j] - fst points.[k]) 
            else 1.0)
        |> List.fold (*) 1.0 
    indizes |> List.sumBy (fun j -> p j * snd points.[j])

示例

这是一个简单的测试 session :

> let points = [|0.0,0.0; 1.0,2.0; 2.0,3.0|];;
val points : (float * float) [] = [|(0.0, 0.0); (1.0, 2.0); (2.0, 3.0)|]

> let f = LagrangeInterpol points;;
val f : (Double -> float)

> f 0.0;;
val it : float = 0.0

> f 1.0;;
val it : float = 2.0

> f 2.0;;
val it : float = 3.0

所以我希望我没有犯任何重大错误。

请注意，我没有在此处进行任何性能优化 - 这应该足以绘制图表或获取中间的一些值。

获取多项式的表示

这有点棘手 - 您可以尝试提出系数的组合公式，或者(像我一样)数学上的懒惰，只需使用足够的运算符实现 Polynom-Type:

type Polynom = 
    Poly of float list with
    override p.ToString () =
        match p with
        | Poly coefs -> 
            System.String.Join (" + ", coefs |> List.mapi (fun i c -> sprintf "%AX^%d" c i))
    static member Const c = Poly [c]
    static member Zero = Polynom.Const 0.0
    static member One = Polynom.Const 1.0
    static member X = Poly [0.0; 1.0]
    static member (+) (Poly cs1, Poly cs2) =
        let m = max (List.length cs1) (List.length cs2)
        List.zip (ofLen m cs1) (ofLen m cs2)
        |> List.map (fun (a,b) -> a+b)
        |> Poly
    static member (-) (Poly cs1, Poly cs2) =
        let m = max (List.length cs1) (List.length cs2)
        List.zip (ofLen m cs1) (ofLen m cs2)
        |> List.map (fun (a,b) -> a-b)
        |> Poly
    static member (*) (f : float, Poly cs2) : Polynom =
        cs2
        |> List.map (fun c -> f * c)
        |> Poly
    static member private shift n (Poly cs) =
        List.replicate n 0.0 @ cs |> Poly
    static member (*) (Poly cs1, p2 : Polynom) : Polynom =
        cs1
        |> List.mapi (fun i c -> Polynom.shift i (c * p2))
        |> List.sum
    static member (/) (Poly cs1, f : float) : Polynom =
        cs1
        |> List.map (fun c -> c / f)
        |> Poly

这里我只是使用一个 float 列表来表示多项式的系数(因此 X^2 + 2X + 3 是 Poly [3.0; 2.0; 1.0]请注意，第 i 系数是位于 X^i 处的系数。

有了这个，我们就可以使用与以前几乎相同的功能:

let getPolynom (points : (float * float)[]) =
    let indizes = [0..points.Length-1]
    let p j = 
        indizes
        |> List.map (fun k -> 
            if k <> j 
            then (Polynom.X - Polynom.Const (fst points.[k]))
                  / (fst points.[j] - fst points.[k]) 
            else Polynom.One)
        |> List.fold (*) Polynom.One
    indizes |> List.sumBy (fun j -> Polynom.Const (snd points.[j]) * p j)      

如您所见，我使用了相同的函数，仅将参数 x 替换为 Polynom.X 并适本地包装了常量。

示例

这里有两个例子(将它们与 Wiki-Page 进行比较，它们应该是正确的):

> LagrangeInterpolation.getPolynom 
    [|(1.0, 1.0); (2.0, 4.0); (3.0, 9.0)|] |> string;;
val it : string = "0.0X^0 + 0.0X^1 + 1.0X^2"

> LagrangeInterpolation.getPolynom 
     [| 1.0,1.0; 2.0,8.0; 3.0,27.0 |] |> string;;
val it : string = "6.0X^0 + -11.0X^1 + 6.0X^2"

带有帮助程序的完整代码

模块内的完整代码是:

module LagrangeInterpolation =

    let private ofLen n cs =
        let l = List.length cs
        if l < n
            then cs @ List.replicate (n-l) 0.0
            else cs

    type Polynom = 
        Poly of float list with
        override p.ToString () =
            match p with
            | Poly coefs -> 
                System.String.Join (" + ", coefs |> List.mapi (fun i c -> sprintf "%AX^%d" c i))
        static member Const c = Poly [c]
        static member Zero = Polynom.Const 0.0
        static member One = Polynom.Const 1.0
        static member X = Poly [0.0; 1.0]
        static member (+) (Poly cs1, Poly cs2) =
            let m = max (List.length cs1) (List.length cs2)
            List.zip (ofLen m cs1) (ofLen m cs2)
            |> List.map (fun (a,b) -> a+b)
            |> Poly
        static member (-) (Poly cs1, Poly cs2) =
            let m = max (List.length cs1) (List.length cs2)
            List.zip (ofLen m cs1) (ofLen m cs2)
            |> List.map (fun (a,b) -> a-b)
            |> Poly
        static member (*) (f : float, Poly cs2) : Polynom =
            cs2
            |> List.map (fun c -> f * c)
            |> Poly
        static member private shift n (Poly cs) =
            List.replicate n 0.0 @ cs |> Poly
        static member (*) (Poly cs1, p2 : Polynom) : Polynom =
            cs1
            |> List.mapi (fun i c -> Polynom.shift i (c * p2))
            |> List.sum
        static member (/) (Poly cs1, f : float) : Polynom =
            cs1
            |> List.map (fun c -> c / f)
            |> Poly

    let getPolynom (points : (float * float)[]) =
        let indizes = [0..points.Length-1]
        let p j = 
            indizes
            |> List.map (fun k -> 
                if k <> j 
                then (Polynom.X - Polynom.Const (fst points.[k]))
                      / (fst points.[j] - fst points.[k]) 
                else Polynom.One)
            |> List.fold (*) Polynom.One
        indizes |> List.sumBy (fun j -> Polynom.Const (snd points.[j]) * p j)      

备注

为了获得更好的输出，您可能应该添加一些简化(例如 Poly [1.0;0.0] -> Poly [1.0])并改进 ToString 方法，但我'我确信你能应付;)