【F#从入门到实战】14. F#实现分部积分法

举报
jackwangcumt 发表于 2021/07/15 10:28:08 2021/07/15
【摘要】 大学高等数学中的微积分是一门重要的基础课程,它广泛应用于数学、物理和工程中。对于求不定积分来说,有一种重要的方法就是分部积分法(integration by parts),它是一种不容易求出积分的情况下,一种很好的尝试方法。本文用F#实现了分部积分法算法。

       欢迎大家来到【F#从入门到实战】,在这里我将分享关于F#编程语言的系列文章,带大家一起去学习和成长,并探索函数编程语言F#这个有趣的世界。所有文章都会结合示例代码和笔者的经验进行讲解,真心想把十余年的IT经验分享给大家,希望对您有所帮助,文章中也定有不足之处,请海涵!本系统文章将从F#基本语法入手,逐步通过自定义类型来实现数学表达式的各种常见解析操作,如对表达式进行求值、化简、展开、求导和求积分等。此系统博文也是了解和实现一个简易的计算机代数系统的基础。

      下面给出【F#从入门到实战】系统专题文章的目录:

【F#从入门到实战】01. F#语言快速入门 
【F#从入门到实战】02. F#数组常见用法   
【F#从入门到实战】03. F#自定义操作符 
【F#从入门到实战】04. F#5.0新特征总结 
【F#从入门到实战】05. F#表达式求值     
【F#从入门到实战】06. F#表达式化简    
【F#从入门到实战】07. F#表达式展开    
【F#从入门到实战】08. F#大整数阶乘    
【F#从入门到实战】09. F#表达式求导   
【F#从入门到实战】10. F#表达式积分   
【F#从入门到实战】11. F#库FParsec入门  
【F#从入门到实战】12. F#库FParsec解析表达式   
【F#从入门到实战】13. F#库FParsec实现求导符号计算   
【F#从入门到实战】14. F#实现分部积分法   

     下面将正式开始本文的介绍: 

     大学高等数学中的微积分是一门重要的基础课程,它广泛应用于数学、物理和工程中。对于求不定积分来说,有一种重要的方法就是分部积分法(integration by parts),它是一种不容易求出积分的情况下,一种很好的尝试方法。关于分部积分法的算法公式如下:

     ∫ u v' dx = uv - ∫ u' v dx

     根据百度百科定义,分部积分法是微积分学中的一类重要的、基本的计算积分的方法。它是由微分的乘法法则和微积分基本定理推导而来的。它的主要原理是将不易直接求结果的积分形式,转化为等价的易求出结果的积分形式的。下面给出一个具体的示例,需要求出如下的不定积分:

2.jpg
     那么分部积分法的基本步骤如下:

3.jpg

       

       那么,如何用F#语言实现这样一种积分算法呢?经过多次尝试,发现一种可行的方法。下面给出核心代码示例。先给出Expr定义:

    type Expr = 
        | CstF of float
        | Sym of string
        | Sin of Expr
        | Cos of Expr
        | Ln of Expr
        | Exp of Expr
        | Pow of Expr * Expr
        | Add of Expr * Expr
        | Sub of Expr * Expr
        | Mul of Expr * Expr
        | Div of Expr * Expr
        | Pi of float
        | Neg of Expr
      
        //~ 一元运算符 - 
        static member (~-) (x:Expr) = 
            Neg(x)
       // 减法
        static member (-) (x:Expr, y:Expr) = 
            Sub(x,y)
        static member (-) (x:Expr, y:float) = 
            Sub(x,CstF y)
        //加法
        static member (+) (x:Expr, y:Expr) = 
            Add(x,y)
        static member (+) (x:float, y:Expr) = 
            Add(CstF x,y)
        static member (+) (x:Expr, y:float) = 
            Add(x,CstF y)
        static member (*) (x:Expr, y:Expr) = 
            Mul(x,y)
         //乘法
        static member (*) (x:float, y:Expr) = 
            Mul(CstF x,y)
        static member (*) (x:Expr, y: float) = 
            Mul(x,CstF y)
        //除法
        static member (/) (x:Expr, y:Expr) = 
            Div(x,y)
        static member (/) (x:Expr, y:float) = 
            Div(x,CstF y)
        static member (/) (x:float, y:Expr) = 
            Div(CstF x,y)
        // ^是特殊符号,这里用 .^ 表示 Power 如 2^3 
        static member (.^) (x:Expr, y:Expr) = 
            Pow(x,y)
        static member (.^) (x:Expr, y:float) = 
            Pow(x,CstF y)
        //e^x
        static member exp (x:Expr) = 
            Exp(x)

     这里的自定义Expr中定义了静态成员方法,是一个自定义操作符,涉及到加减乘除,且支持不同的重载(Overload) 。由于分部积分涉及到求导和基本的积分规则,下面给出定义:

 
    let rec diff (expr : Expr , v ) = 
        match expr with
        | CstF x  -> CstF 0.
        | Sym x when x =  v -> CstF 1.
        | Sym x when x <>  v -> CstF 0.
        | Sin(x) when x = Sym v -> Cos(x)
        | Sin(x)  -> Cos(x) * diff(x,v)
        | Cos(x) when x = Sym v -> -Sin(x)
        | Cos(x)  -> -Sin(x) * diff(x,v)
        | Pow(a,x) when x = Sym v -> Ln(a) * Pow(a,x)
        | Ln(x) when x = Sym v -> 1. / x
        | Pow(x,y) when x = Sym v -> y * Pow(x, y - 1.)
        | Exp(x) when x = Sym v -> Exp(x)
        | Exp(x)  -> Exp(x) * diff(x,v)
        | Add(x,y) -> diff(x,v) + diff(y,v) //(f(x)+g(x))'=f'(x)+g'(x)
        | Sub(x,y) -> diff(x,v) - diff(y,v) //(f(x)-g(x))'=f'(x)-g'(x)
        | Mul(CstF n,f) ->  n * diff(f,v)
        | Mul(f,g) ->  diff(f,v)* g - f * diff(g,v) //(f(x)g(x))'=f'(x)g(x)+f(x)g'(x)
        | Div(f,g) -> (diff(g,v) * f -  diff(f,v) * g)/ (f .^ 2.) //(g(x)/f(x))'=(g'(x)f(x)-f'(x)g(x))/(f(x))^2
        | e -> e



    let rec simp expr = 
        match expr with
        | CstF f  -> CstF f 
        | Sym x  -> Sym x
        | Neg(Neg(e)) -> simp e
        | Add(CstF f ,CstF f2) -> CstF (f+f2)
        | Sub(CstF f ,CstF f2) -> CstF (f-f2)
        | Mul(CstF f ,CstF f2) -> CstF (f * f2)
        | Mul(CstF f ,Mul(CstF f2,e)) -> simp (Mul(CstF (f * f2),simp e))
        | Div(CstF f ,CstF f2) -> CstF (f / f2)
        | Div(Div(e, CstF f), CstF f2) -> Div(simp e,CstF (f * f2))
        | Sub(e1 ,Neg(e2)) -> Add(e1,e2)
        | Sub(x,y) -> Sub(simp x,simp y)
        | Add(CstF 0. ,e) -> simp e
        | Add(e,CstF 0. ) -> simp e
        | Add(x,y) -> Add(simp x,simp y)
        | Mul(CstF 1.,e) -> simp e
        | Mul(CstF 0.,e) -> CstF 0.
        | Mul(e,CstF 0.) -> CstF 0.
        | Mul(x,y) -> Mul(simp x,simp y)
        | e -> e

    let rec int2 (expr : Expr , v ) = 
        match expr with
        | CstF f  -> Sym v 
        | Sym v -> 0.5 * Pow(Sym v, CstF 2.)
        | Sin(x) when x = Sym v -> -Cos(x)
        | Cos(x) when x = Sym v -> Sin(x)
        | Cos (Div (x, CstF n)) when x = Sym v -> n * Sin (Div (x, CstF n))
        | Sin (Div (x, CstF n)) when x = Sym v -> - n * Cos (Div (x, CstF n))
        | Exp(x) when x = Sym v -> Exp(x)
        | Exp(Mul(CstF n,x)) when x = Sym v -> Exp(Mul(CstF n,x)) / n
        | Exp(Mul(x,CstF n)) when x = Sym v -> Exp(Mul(CstF n,x)) / n
        | Add(x,y) -> int2(x,v) + int2(y,v)
        | Sub(x,y) -> int2(x,v) - int2(y,v)
        | Mul(CstF 1.,e) -> int2(e,v)
        | Mul(CstF n,e) -> n * int2(e,v)
        | Pow(x,y) when x = Sym v ->  Pow(x, y + 1.) / ( y + 1.)
        | Div(e,CstF 1.) -> int2(e,v)
        | Div(e,CstF n) -> int2(e,v) / n
        | e -> e

    let rec int (expr : Expr , v ) = 
        match expr with
        | CstF f  -> Sym v 
        | Sym v -> 0.5 * Pow(Sym v, CstF 2.)
        | Sin(x) when x = Sym v -> -Cos(x)
        | Cos(x) when x = Sym v -> Sin(x)
        | Cos (Div (x, CstF n)) when x = Sym v -> n * Sin (Div (x, CstF n))
        | Sin (Div (x, CstF n)) when x = Sym v -> - n * Cos (Div (x, CstF n))
        | Exp(x) when x = Sym v -> Exp(x)
        | Exp(Mul(CstF n,x)) when x = Sym v -> Exp(Mul(CstF n,x)) / n
        | Exp(Mul(x,CstF n)) when x = Sym v -> Exp(Mul(CstF n,x)) / n
        | Add(x,y) -> int(x,v) + int(y,v)
        | Sub(x,y) -> int(x,v) - int(y,v)
        | Pow(x,y) when x = Sym v ->  Pow(x, y + 1.) / ( y + 1.)
        | Mul(CstF 1.,e) -> int(e,v) //防止无限循环
        | Mul(CstF n,e) -> n * int(e,v) //防止无限循环
        | Div(e,CstF 1.) -> int2(e,v)
        | Div(e,CstF n) -> int2(e,v) / n
        | Mul(u,t) -> 
            printfn "%O" u
            printfn "%O" t
            let tf = simp (int2(t,v))
            printfn "%O" tf
            let du = simp (diff(u,v))
            let z = u * tf - int2(du*tf ,v) //分部积分法
            z
        | e -> e

   其中的分部积分规则匹配一个uv的乘法,通过反求积分后,对算法进行改造,即 z = u * int(t,x) - int(diff(u,x)*int(t,x) ,x) 。这里省略了打印的函数定义。另外,需要注意,代码中定义了一个两个非常相似的积分函数,即int和int2。其中差异就是int2中未有分部积分定义,这样可以防止无限循环。

    let x = Sym "x"
    let y =  (CstF 6.) * x // 6*x
    let f = x * Exp(y)
    let z = int(f,"x")
    printfn "%s" (simp z |> printExpr)

    在主程序中,需要注意的就是6x不能直接写成 6.*x,而需要写成 (CstF 6.) * x ,这个好像Expr里面定义的乘法重载未生效,只生效了Expr * Expr这种形式的,具体原因不明,可能是内部的自定义操作符,不是全局定义的。 执行后结果如下,这与开头的示例一致:

((x*(exp((6*x))/6))-(exp((6*x))/36))

    同理,再给出一个示例:

    let f = x * Cos(x)
    let z = int(f,"x")
    printfn "%s" (simp z |> printExpr)

   输出结果如下:

(x*sin(x))+cos(x)

   在给出一个相对复杂一点的案例:

5.jpg

    其解题步骤为:

6.jpg

   

    let t = Sym "t"
    let y = (CstF 3. * t + CstF 5. )* Cos( t / CstF 4. )
    let z = int(y,"t")
    printfn "%s" (simp z |> printExpr)

    执行代码,结果如下:

((((3*t)+5)*(4*sin((t/4))))-(-48*cos((t/4))))

    化简后为:

12*sin(t/4)*t+20*sin(t/4)+48*cos(t/4)

     这与上述示例结果一致。注意,如果其他的积分求的不对,则需要补充特定规则的积分和求导规则,让其分部积分可以正确处理。

【版权声明】本文为华为云社区用户原创内容,转载时必须标注文章的来源(华为云社区)、文章链接、文章作者等基本信息, 否则作者和本社区有权追究责任。如果您发现本社区中有涉嫌抄袭的内容,欢迎发送邮件进行举报,并提供相关证据,一经查实,本社区将立刻删除涉嫌侵权内容,举报邮箱: cloudbbs@huaweicloud.com
  • 点赞
  • 收藏
  • 关注作者

评论(0

0/1000
抱歉,系统识别当前为高风险访问,暂不支持该操作

全部回复

上滑加载中

设置昵称

在此一键设置昵称,即可参与社区互动!

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。