【F#从入门到实战】15. F#实现换元积分法

举报
jackwangcumt 发表于 2021/07/23 08:13:53 2021/07/23
【摘要】 大学高等数学中的微积分是一门重要的基础课程,它广泛应用于数学、物理和工程中。对于求不定积分来说,除了上一篇博文介绍的分部积分法(integration by parts)外,还有一种常见的积分方法,即积分换元法。

       欢迎大家来到【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#实现分部积分法   

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

1 换元法概述

     大学高等数学中的微积分是一门重要的基础课程,它广泛应用于数学、物理和工程中。对于求不定积分来说,除了上一篇博文介绍的分部积分法(integration by parts)外,还有一种常见的积分方法,即积分换元法。下面给出一个实际的数学换元示例:

1.jpg
      如果需要求上述表达式的不定积分,则直接进行求解有点难度,此时可以尝试使用换元法进行求解,这个表达式的分母比较复杂,将分母作为切入口进行换元,具体的换元过程如下:

2.jpg
      此示例来自https://tutorial.math.lamar.edu/Classes/CalcI/SubstitutionRuleIndefinite.aspx 。此网站还有其他的数学知识,非常好,感兴趣的可以阅读一下。

2 F#实现换元法积分

      换元法的过程比较灵活,有时候需要通过凑的方式才能顺利求出积分。因此,完成一个功能完备的换元积分程序是不容易的。本文只是针对特定几个场景进行换元法尝试,主要是给出一种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
        | Sqrt of Expr
        | D of Expr

       其次,给出一个求微分的函数,代码如下:

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

      再次,换元法则需要将表达式当中的特定对象进行替换,给出替换函数,代码如下:

    let rec substitute (expr : Expr , x : Expr , u : Expr ) = 
        match expr with
        //| x  -> u //有问题
        | e when e = x -> u
        | Cos(e) when e = x -> Cos(u)
        | Sin(e) when e = x-> Sin(u)
        | Exp(e) when e = x-> Exp(u)
        | Pow(e,n) when e = x -> Pow(u,n)
        | Add(e1,e2) -> Add(substitute(e1,x,u),substitute(e2,x,u))
        | Sub(e1,e2) -> Sub(substitute(e1,x,u),substitute(e2,x,u))
        | Mul(e1,e2) -> Mul(substitute(e1,x,u),substitute(e2,x,u))
        | Div(e1,e2) -> Div(substitute(e1,x,u),substitute(e2,x,u))
        | Neg(e1) -> Neg(substitute(e1,x,u))
        | Ln(e1) -> Ln(substitute(e1,x,u))
        | e -> e

       这里需要注意一下,F#中的模式匹配当中,匹配的项不能直接引用函数参数,否则可能会出现不匹配的问题。如替换函数中的参数有一个x参数,那么在匹配规则中如果直接用x,如Sin(x),那么则可能出现不匹配的情况,而是用Sin(e) when e = x 。第二个就是匹配规则中,不能用同一个变量,如Add(Sin(e),e),而需要用Add(Sin(e1),e2) 。
      再次,再给出一个dsolve辅助函数,这里主要就是求解类似于  y*dy = du 反求dy = du / y 的过程,这个函数非常初级,如果需要完善,则需要构建更加多的,覆盖面广的规则:

    let rec dsolve (du : Expr , dt : Expr ) = 
        match du with
        | Mul(e1,dy) when dt = dy -> Div(CstF 1.,e1)
        | Mul(e1,e2) -> Mul(Div(CstF 1.,e1), dsolve(e2,dt))
        | e -> e

     最后,给出换元法的核心代码:

   //换元法
    let rec int_subs (expr : Expr , v ) = 
        match expr with
        | Add(e1,e2) -> Add(int_subs(e1,v),int_subs(e2,v))
        | Sub(e1,e2) -> Sub(int_subs(e1,v),int_subs(e2,v))
        | Cos(e)  -> 
            printfn "%O" e
            match e with
            | Sqrt(Sym v) -> 
                //printfn "%O" v
                let u = Sym "u"
                let x = Pow(u,CstF 2.)
                let dx = simp2 (diff(x,"u"))
                //printfn "%O" dx
                let st = simp2 (substitute(expr,e,u))
                let st = int_subs (dx * st , "u")
                let st = simp2 st
                let stf = substitute(st,u,e)
                stf
            | e -> e   
        | Exp(e) -> 
            printfn "%O" e
            match e with
            | Sqrt(Sym v) -> 
                //printfn "%O" v
                let u = Sym "u"
                let x = Pow(u,CstF 2.)
                let dx = simp2 (diff(x,"u"))
                //printfn "%O" dx
                let st = simp2 (substitute(expr,e,u))
                let st = int_subs (dx * st , "u")
                let st = simp2 st
                let stf = substitute(st,u,e)
                stf
            | e -> e   
        | Div(e1,e2) -> 
            //printfn "%O" e1
            //printfn "%O" e2
            let u = e2
            let du = simp2 (d(u,Sym v)) // du = u'dv
            printfn "%O" du
            let dv = dsolve(du,d(Sym v,Sym v))
            printfn "%O" (simp2 dv)
            let u = Sym "u"
            let st =simp2 (substitute (expr ,e2, u))
            printfn "%O" (simp2 st)
            let z =simp2 (st  * dv)
            printfn "%O" (simp2 z)
            let f = simp2(int2(z,"u"))
            printfn "%O" (f)
            let fs =simp2 (substitute (f ,u,e2))
            fs
        | e -> e

    注意:这里的换元法只处理了少数几个场景。另外,关于化简函数simp2 以及 打印函数这里省略。感兴趣的可以查看之前的相关博文。

3 代码测试

      下面给出测试代码,来验证换元法是否可以按预期进行运行,其代码如下:

    let y = Sym("y")
    let z = (CstF 3. * y) / (CstF 5. * (y .^2.) + CstF 4. )
    let z = int_subs (z,"y")
    printfn "%O" z
    printfn "%s" (simp2 z |> printExpr)
    ///////////////////////////////////////////////////////////
    let x = Sym("x")
    let y = Cos(Sqrt(x))
    let z = int_subs (y,"x")
    //printfn "%O" z
    printfn "%s" (simp2 z |> printExpr)
    ///////////////////////////////////////////////////////////
    let y = Exp(Sqrt(x))
    let z = int_subs (y,"x")
    //printfn "%O" z
    printfn "%s" (simp2 z |> printExpr)

     运行结果如下:

((3/10)*ln(((5*(y^2))+4)))
((2*(sqrt(x)*sin(sqrt(x))))+(2*cos(sqrt(x))))
((2*(sqrt(x)*exp(sqrt(x))))-(2*exp(sqrt(x))))

    对括号进行化简后如下:

(3/10)*ln(5*y^2+4)
2*sin(sqrt(x))*sqrt(x)+2*cos(sqrt(x)) 
2*sqrt(x)*e^sqrt(x)-2*e^sqrt(x)

   即对应如下公式:

6.jpg

6.jpg

7.jpg

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

评论(0

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

全部回复

上滑加载中

设置昵称

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

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

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

举报
请填写举报理由
0/200