Аналитическое вычисление производной функции на языке Scala11.09.2016 02:03

Введение


Данный алгоритм реализован на языке 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 — переменная, по которой происходит дифференцирование.

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

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), возвращаемое значение — функция-производная 08d2849d88a2cb274361121b1a61ef7c.gif в том же формате.

Примечание 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 — модуль, поскольку его нет в таблице.

4cfe0ca940ca2212d4fdfeba15bbcbcb.gif

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


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


  • Математический портал «mathprofi.ru» mathprofi.ru/slozhnye_proizvodnye_logarifmicheskaja_proizvodnaja.html
  • Odersky M. «Programming in Scala (second edition)»
  • Хостманн К. «Scala для нетерпеливых»

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

© Habrahabr.ru