Аналитическое вычисление производной функции на языке Scala

    Введение


    Данный алгоритм реализован на языке Scala, характерной особенностью которого является использование case-классов, так удачно подходящих для написания алгоритма дифференцирования. В этой статье планируется описать лишь часть программы, содержащей алгоритм нахождения производной, поскольку разработка парсера для математических выражений это другая большая тема,
    заслуживающая отдельной статьи


    Подготовка


    Сначала опишем структуру данных, в которой будет храниться исходная математическая функция. Опишем трейт MathAST:

    sealed trait MathAST
    

    И его наследников:

    sealed trait SingleToken extends MathAST{
    	val a: MathAST
    }
    sealed trait DoubleToken extends MathAST{
    	val a: MathAST
    	val b: MathAST
    }
    sealed trait LeafToken extends MathAST
    

    SingleToken будет реализовывать case-классы унарных операторов, таких как sin(a), -a, ln(a) и т.д.
    DoubleToken описывает бинарные операторы: a+b, a^b и т.д.
    LeafToken — «листья» дерева, т.е. константы, переменные и зарезервированные имена констант (число Пи и экспонента).

    Опишем классы/объекты операторов и токенов:

    case object Pi extends LeafToken
    case object Exponenta extends LeafToken
    
    case class Sin(override val a: MathAST) extends SingleToken
    case class Cos(override val a: MathAST) extends SingleToken
    …
    case class Mul(override val a: MathAST, override val b: MathAST) extends DoubleToken
    case class Add(override val a: MathAST, override val b: MathAST) extends DoubleToken
    …
    case class Differentiate(f: MathAST, dx: Variable) extends MathAST
    
    case class Variable(name: String) extends LeafToken
    case class Constant(value: BigDecimal) extends LeafToken

    Обратите внимание на класс Differentiate, он имеет особую сигнатуру: f – исходная функция, dx – переменная, по которой происходит дифференцирование.

    Теперь есть все, чтобы представить математическую функцию в виде дерева вычислений, для примера возьмем функцию: , которая примет вид:

    Mul(Constant(BigDecimal(2)), Pow(x, Constant(BigDecimal(2)))

    Конечно, чтобы получить дерево-выражение из обычной строки, введенной пользователем, нужно написать парсер, но, как было упомянуто выше, это уже другая тема. Скажу лишь что в программе используется самодельный парсер, наследующий трейт Parsers из пакета scala.util.parsing.combinator.


    Алгоритм нахождения производной


    В основе которого лежат правила дифференцирования и таблица производных.

    Опишем рекурсивную функцию, которая и будет преобразовывать исходную математическую функцию в результирующую функцию-производную:

    def differentiate(f: MathAST)(implicit dx: String): MathAST

    Аргумент dx, содержащий имя переменной (по которой происходит дифференцирование) помечен как неявный (implicit), это позволит не передавать ее в рекурсивные вызовы, пусть этим занимается компилятор.

    На вход рекурсивной функции подается выражение — исходная функция f(x) (в формате MathAST), возвращаемое значение — функция-производная в том же формате.

    Примечание 1: Выражение может быть бинарным, унарным или токеном.
    Примечание 2: Оператором может быть один из: «+», «-», «^», «*», «/», «abs», «sin», «cos», «tg», «ctg», «ln», «arcsin», «arccos», «arctg», «arcctg», «(», «)».
    Примечание 3: Входные и выходные данные представлены в формате MathAST — дерево-выражение.


    Общий алгоритм


    В общем виде алгоритм слишком абстрактный, поэтому дальше разберем его подробнее.


    1. Рекурсивная функция получает на вход данные и используя сопоставление с образцом (pattern-matching) выполняет необходимые действия, в зависимости от типа данных.

    2. Функция высчитывает производную для входного выражения и возвращает выражение-результат. Может получиться, что аргументы a и/или b оказались не константой и не переменной, а сложной функцией u(x),
      тогда понадобится рекурсивно посчитать еще и производную u’(x), т.е. вернуть [ differentiate(u(x)) ] — перейти к шагу 1 с новыми данными — u(x).

    3. Если данные не корректны вернуть сообщение об ошибке.

    Детали принципа работы и связь с математическими абстракциями


    Функция приняла на вход данные — выражение-функцию, которую следует обработать в соответствии с правилами дифференцирования


    Если бинарное выражение


    Бинарные выражения помечены трейтом DoubleToken. Оператор проверяется на соответствие одному из доступных операторов («+», «-», «*», «/», «^»):

    case Add(a, b) => Add(differentiate(a), differentiate(b))
    case Sub(a, b) => Sub(differentiate(a), differentiate(b))
    …
    


    1. Если оператор «+»: вернуть [ differentiate(a) + differentiate(b) ].

    2. Если оператор «-»: вернуть [ differentiate(a) — differentiate(b) ].

    3. Если оператор «*»: Умножение представляет из себя более сложный случай, операнды a и b могут быть константами или переменными (всего 4 комбинации: u(x)*c, u(x)*v(x), c*c, c*u(x)).
      Функция анализирует какой из 4 вариантов попался и возвращает выражение используя правило дифференцирования № 1, № 3, и №5,
      если один из операндов – сложная функция. Например: если a = u(x), а b = v(x), то вернуть [ differentiate(a) * b + a * differentiate(b)) ].
      Приватный метод isDependsOnVar проверяет, зависит ли подвыражение от переменной, по которой производится дифференцирование.

    4. Если оператор «/»: Действия аналогичны, как и с оператором «*».

    5. Если оператор возведения в степень «^»: Здесь имеется так же 4 варианта выражений: u(x)^c, u(x)^v(x), c^c, c^u(x).
      Третий случай самый нетривиальный: u(x)^v(x), чтобы найти производную от такого выражения нужно использовать логарифм и его свойства (описание процедуры логарифмирования выходит за рамки этой статьи, поэтому сразу приведем готовую формулу – № 6).

      В остальных трех случаях используются стандартные табличные формулы.

    В качестве примера приведем обработку операции умножения:

    
    case Mul(a, b) => {
    	val u = isDependsOnVar(a)
    	val v = isDependsOnVar(b)
    	if (u && !v) {
    		Mul(differentiate(a), b)
    	} else if (!u && v){    // c^u(x), c=const
    		Mul(a, differentiate(b))
    	}else if (!u && !v){    // c^c, c=const
    		Constant(BigDecimal(0))
    	}else Add(Mul(differentiate(a), b), Mul(a, differentiate(b)))
    }


    Если унарное выражение


    Классы SingleToken обрабатываются следующим образом:

    
    case e: SingleToken =>
    val d = e match {
    	case Sin(x) => Cos(x)
    	case Cos(x) => Usub(Sin(x))
    	case Tg(x) => Div(Constant(BigDecimal(1)), Pow(Cos(x), Constant(BigDecimal(2))))
    	…
    }
    if (isLeaf(e.a)) d else Mul(d, differentiate(e.a))

    Оператор проверяется на соответствие одному из доступных операторов («sin», «-», «cos», …)
    Для примера, оператор «sin»: вернуть [ cos(a) ], если a = x, если же a — сложная функция u(x), то вернуть [ cos(a) * differentiate(a) ].

    С остальными операторами происходят аналогичные действия, используя правило дифференцирования сложной функции и табличные правила взятия производной.

    Отдельно следует рассмотреть оператор abs — модуль, поскольку его нет в таблице.


    Приватный метод isLeaf выясняет сложная функция или нет, в первом случае нужно вернуть текущий результат умноженный на производную вложенной функции, а во втором просто вернуть текущий результат.


    Если один токен


    Речь идет о переменной или константе

    case Variable(a) => if (a == dx) Constant(BigDecimal(1)) else Constant(BigDecimal(0))
    case Constant(a) => Constant(BigDecimal(0))
    case Pi | Exponenta => Constant(BigDecimal(0))
    case _ => throw new AbstractEvaluateException("Differentiate: Wrong input data")


    1. Введены некорректные данные, вывести сообщение об ошибке и завершить работу.
    2. Если переменная (по которой осуществляется дифференцирование, например x), вернуть [ 1 ].
    3. Если константа, вернуть [ 0 ].

    Напоследок добавим строку:

    case Differentiate(_f, _dx) => differentiate(_f)(_dx.name)

    Это на случай, если внутри функции есть вложенная производная, т.е. теперь можно считать производные высших порядков.

    Приведу весь код дифференцирования, для лучшего понимания:

    object Derivative {
      def apply(e: MathAST)(implicit dx: String): MathAST = differentiate(e)
    
      def differentiate(f: MathAST)(implicit dx: String): MathAST = f match {
        case Differentiate(_f, _dx) => differentiate(_f)(_dx.name)
    
        case Add(a, b) => Add(differentiate(a), differentiate(b))
        case Sub(a, b) => Sub(differentiate(a), differentiate(b))
        case Div(a, b) => Div(Sub(Mul(differentiate(a), b), Mul(a, differentiate(b))), Pow(b, Constant(BigDecimal(2))))
    
        case Pow(a, b) => {
          val u = isDependsOnVar(a)
          val v = isDependsOnVar(b)
          if (u && !v) Mul(Mul(Pow(a, Sub(b, Constant(BigDecimal(1)))), differentiate(a)), b) // u(x)^c, c=const
          else if (!u && v){  // c^u(x), c=const
            Mul(Mul(Pow(a, b), differentiate(b)), Ln(a))
          }else if(!u && !v){
            Constant(BigDecimal(0))
          }else Mul(Pow(a, Sub(b, Constant(BigDecimal(1)))), Add(Mul(b, differentiate(a)), Mul(Mul(a, Ln(a)), differentiate(b))))
        }
        case Mul(a, b) => {
          val u = isDependsOnVar(a)
          val v = isDependsOnVar(b)
          if (u && !v) {
            Mul(differentiate(a), b)
          } else if (!u && v){  // c^u(x), c=const
            Mul(a, differentiate(b))
          }else if (!u && !v){// c^c, c=const
            Constant(BigDecimal(0))
          }else Add(Mul(differentiate(a), b), Mul(a, differentiate(b)))
        }
    
        case e: SingleToken =>
          val d = e match {
            case Sin(x) => Cos(x)
            case Cos(x) => Usub(Sin(x))
            case Tg(x) => Div(Constant(BigDecimal(1)), Pow(Cos(x), Constant(BigDecimal(2))))
            case Ctg(x) => Usub(Div(Constant(BigDecimal(1)), Pow(Sin(x), Constant(BigDecimal(2)))))
            case Abs(x) => Div(x, Abs(x))
            case Ln(x) => Div(Constant(BigDecimal(1)), x)
            case Sqrt(x) => Div(Constant(BigDecimal(1)), Mul(Constant(BigDecimal(2)), Sqrt(x)))
            case Usub(x) => Usub(differentiate(x))
            case Arcsin(x) => Div(Constant(BigDecimal(1)), Sqrt(Sub(Constant(BigDecimal(1)), Pow(x, Constant(BigDecimal(2))))))
            case Arccos(x) => Usub(Div(Constant(BigDecimal(1)), Sqrt(Sub(Constant(BigDecimal(1)), Pow(x, Constant(BigDecimal(2)))))))
            case Arctg(x) => Div(Constant(BigDecimal(1)), Sub(Constant(BigDecimal(1)), Pow(x, Constant(BigDecimal(2)))))
            case Arcctg(x) => Usub(Div(Constant(BigDecimal(1)), Sub(Constant(BigDecimal(1)), Pow(x, Constant(BigDecimal(2))))))
            case _ => throw new AbstractEvaluateException("Differentiate: Unknown unary operator")
          }
          if (isLeaf(e.a)) d else Mul(d, differentiate(e.a))
    
        case Variable(a) => if (a == dx) Constant(BigDecimal(1)) else Constant(BigDecimal(0))
        case Constant(a) => Constant(BigDecimal(0))
        case Pi | Exponenta => Constant(BigDecimal(0))
        case _ => throw new AbstractEvaluateException("Differentiate: Wrong input data")
      }
    
      private def isLeaf(e: MathAST): Boolean = e match {
        case Variable(_) | Constant(_) => true
        case  Pi | Exponenta => true
        case _ => false
      }
    
      private def isDependsOnVar(tree: MathAST)(implicit dx: String): Boolean = tree match{
        case e: DoubleToken => (e.a match {
          case Variable(name) => if(name == dx) true else false
          case _ => isDependsOnVar(e.a)
        })||(e.b match {
          case Variable(name) => if(name == dx) true else false
          case _ => isDependsOnVar(e.b)
        })
        case e: SingleToken => isDependsOnVar(e.a)
        case Variable(name) => if(name == dx) true else false
        case _ => false
      }
    }


    Заключение


    Весь код исходников можно скачать на github'е, протестировать программу онлайн можно на сайте Калькулятор производных онлайн, приложение выполнено в виде REST сервиса и дополнено модулями упрощения выражений.


    Список литературы


    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More
    Ads

    Comments 13

      0

      А почему не на макросах? case классы — это же элементарно и совершенно очевидно.
      Некоторые даже на С++ шаблонах дифференцирование делают.

        +3
        Символьное вычисление — очень неэффективный способ автоматического дифференцирования, на практике применим разве что для символьного решения школьных заданий. (Не критика, т.к. в заголовке сказано, что способ аналитический, но совет для тех, кому нужны автоматические производные для реальных вычислительных задач). В реальности используются совсем другие способы, гораздо более эффективные. Хороший пример реализации метода с дуальным числом на C++ — в гугловской библиотеке Ceres Solver.
          0
          Спасибо за ссылки, не знал о таком
          –1
          Интересно. Но на x^x онлайн-калькулятор даёт неверный результат.
            0
            d(x^x)/d(x) выдает (x^(x-1))*(x+ln(x)) = (x^x)*(1+ln(x)), все верно
            0
            почему бы функциям самим не провайдить способ своего дифференцирования?
              0
              подумаю над этим, спасибо
              +1
              Вообще дифференцирование произведения (а также частного или степени) не зависит о того константы указанные функции или нет. Формула (uv)' = u'v + v'u верна всегда. Поэтому можно было отдельные случаи не расписывать. Просто умножение на 0 надо корректно обработать.
                0

                Как вы сказали, формула (uv)' = u'v + v'u всегда корректна, а после вычисления производной нужно просто ввести дополнительную фазу — упрощение выражения. Ее также можно реализовать с помощью описанного в статье подхода (применять такие правила, как a + 0 = a, a 1 = a, a b = b * a и т.д.). Я тоже таким занимался, только на C#. Почитать можно здесь: Упрощение выражений.

                0
                Можно было не писать парсер а использовать перезагрузку операторов которая в scala из коробки. Ну где-то как это реализовано в scipy.
                  0
                  А если уже откомпилированной программу надо взять выражение со стороны — использовать аналог eval?
                    0
                    А это уже другая задача)
                    Я б написал библиотеку с возможностью юзать ее в repl-е.
                    Дальше такую вещь можно много куда еще интегрировать.
                    Но если вам нравится парсить — то это таки ваше право)
                  0
                  А почему скала? То есть я вижу что это скала ноге же все плюшки — фуксии высших порядков, карирование, отложены вычисления.

                  Only users with full accounts can post comments. Log in, please.