10.09.2016 21:43
mathsatan 13 5500 Источник

Введение


Данный алгоритм реализован на языке 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 сервиса и дополнено модулями упрощения выражений.


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


Поделиться с друзьями
-->

Комментарии (13)


  1. Akon32
    11.09.2016 09:05
    #9800066

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


  1. vagran
    11.09.2016 10:34
    #9800106
    +3

    Символьное вычисление — очень неэффективный способ автоматического дифференцирования, на практике применим разве что для символьного решения школьных заданий. (Не критика, т.к. в заголовке сказано, что способ аналитический, но совет для тех, кому нужны автоматические производные для реальных вычислительных задач). В реальности используются совсем другие способы, гораздо более эффективные. Хороший пример реализации метода с дуальным числом на C++ — в гугловской библиотеке Ceres Solver.


    1. mathsatan
      11.09.2016 12:48
      #9800230

      Спасибо за ссылки, не знал о таком


  1. ik62
    11.09.2016 11:42
    #9800150
    -1

    Интересно. Но на x^x онлайн-калькулятор даёт неверный результат.


    1. mathsatan
      11.09.2016 12:41
      #9800218

      d(x^x)/d(x) выдает (x^(x-1))*(x+ln(x)) = (x^x)*(1+ln(x)), все верно


  1. fogone
    11.09.2016 12:31
    #9800206

    почему бы функциям самим не провайдить способ своего дифференцирования?


    1. mathsatan
      11.09.2016 12:49
      #9800234

      подумаю над этим, спасибо


  1. bay73
    11.09.2016 13:27
    #9800262
    +1

    Вообще дифференцирование произведения (а также частного или степени) не зависит о того константы указанные функции или нет. Формула (uv)' = u'v + v'u верна всегда. Поэтому можно было отдельные случаи не расписывать. Просто умножение на 0 надо корректно обработать.


    1. KvanTTT
      11.09.2016 13:54
      #9800288

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


  1. AStek
    11.09.2016 14:52
    #9800358

    Можно было не писать парсер а использовать перезагрузку операторов которая в scala из коробки. Ну где-то как это реализовано в scipy.


    1. potan
      12.09.2016 14:07
      #9801740

      А если уже откомпилированной программу надо взять выражение со стороны — использовать аналог eval?


      1. AStek
        12.09.2016 17:19
        #9802112

        А это уже другая задача)
        Я б написал библиотеку с возможностью юзать ее в repl-е.
        Дальше такую вещь можно много куда еще интегрировать.
        Но если вам нравится парсить — то это таки ваше право)


  1. kafalk
    14.09.2016 11:56
    #9805564

    А почему скала? То есть я вижу что это скала ноге же все плюшки — фуксии высших порядков, карирование, отложены вычисления.