Символьное дифференцирование

ОГЛАВЛЕНИЕ

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

•    Скачать исходники - 36.24 Кб (v0.90071)
•    Скачать демонстрационный проект - 170.4 Кб (v0.90071)

Введение

Приложение, как на рисунке выше, содержит простую сетку для рисования входной кривой и взятия ее производной на одном и том же виде. Детальное описание математики выходит за рамки статьи. В следующем разделе дано только определение производной. Если вам не интересна математика, пропустите его.

Определение дифференцирования

Производную функции можно трактовать как функцию, значение которой в x равняется градиенту касательной к графику y=f(x) в x, или же как функцию, описывающую мгновенную скорость изменения y относительно x в точке x. Производную f относительно x также можно обозначить через

d(f(x))/dx,
or if y = f(x), dy/dx

Формулы дифференцирования

Выражение

Производная

Замечания

 

u+v

du/dx+dv/dx

 
 

u-v

du/dx-dv/dx

 
 

u/v

(v*du/dx-u*dv/dx)/v^2

 
 

u*v

u*dv/dx+v*du/dx

 
 

c*u

c*du/dx

c - константа

 

u^v

v*u^(v-1)*du/dx+u^v*ln(u)*dv/dx

 
 

u^n

n*u^(n-1)*du/dx

n- вещественный

 

c^u

c^u*ln(c)*du/dx

c - константа

 

e^u

e^u*du/dx

e = 2.7182818284590452353602874713527

1

sin(u)

cos(u)*du/dx

 

2

cos(u)

-sin(u)*du/dx

 

3

tan(u)

sec(u)^2*du/dx

 

4

sec(u)

sec(u)*tan(u)*du/dx

 

5

cosec(u)

-cosec(u)*cot(u)*du/dx

 

6

cot(u)

-cosec(u)^2*du/dx

 

7

sinh(u)

cosh(u)*du/dx

 

8

cosh(u)

sinh(u)*du/dx

 

9

tanh(u)

sech(u)^2*du/dx

 

10

sech(u)

sech(u)*tanh(u)*du/dx

 

11

cosech(u)

cosech(u)*coth(u)*du/dx

 

12

coth(u)

-cosech(u)^2*du/dx

 

13

asin(u)

1/sqrt(1-u^2)*du/dx

 

14

acos(u)

-1/sqrt(1-u^2)*du/dx

 

15

atan(u)

1/(1+u^2)*du/dx

 

16

asec(u)

1/(|u|*sqrt(u^2-1))*du/dx

|u| - abs(u)

17

acosec(u)

-1/(|u|*sqrt(u^2-1))*du/dx

|u| - abs(u)

18

acot(u)

-1/(1+u^2)*du/dx

 

19

asinh(u)

1/sqrt(u^2+1)*du/dx

 

20

acosh(u)

1/sqrt(u^2-1)*du/dx

 

21

atanh(u)

1/(1-u^2)*du/dx

 

22

asech(u)

-1/(u*sqrt(1-u^2))*du/dx

 

23

acosech(u)

-1/(u*sqrt(1+u^2))*du/dx

 

24

acoth(u)

1/(1-u^2)*du/dx

 

25

sqrt(u)

1/(2*sqrt(u))*du/dx

 

26

log10(u)

1/(u*ln(10))*du/dx

 

27

log(u)

1/u*du/dx

 

28

ln(u)

1/u*du/dx

 

29

sign(u)

0

 

30

abs(u)

u/|u|*du/dx

|u| -abs(u)

Математическое выражение

Математическое выражение содержит числа, операторы, функции и переменные. Операторы включают в себя +, -, *, / и ^. Функции подобны sin(x), а переменные подобны x. Задача – определить каждый элемент и применить правильное правило дифференцирования к этому элементу или группе элементов.

Примеры выражений:

sin(2*x)/x
sin(45+cos(2)/tan(x)
sign(cos(x)
32*9-8/2

Шаги дифференцирования

Разбор выражения для заполнения набора

Разбор выражения – процесс сканирования выражения для отыскания его элементов (чисел, операторов, функций и переменных). Каждый оператор должен иметь два операнда; например, 2*12. Оператор "*" имеет два операнда: 2 и 12. Каждая функция имеет аргумент. Например, для sin(x) функцией является sin, а аргументом - x. Надо особым образом сканировать выражение, чтобы определить каждый оператор и операнды, а также каждый аргумент функции, чтобы упростить расчет или дифференцирование.

Выполнение выражения зависит от приоритета операторов. Если построить такой набор операторов и операндов, чтобы за каждым оператором следовали его два операнда, то 2*12 построит следующий набор: "*, "2", "12". Итак, функция расчета проверяет набор. Если элемент является оператором, то функция берет два его оператора и производит расчет для них по рекурсивной формуле, чтобы вычислить все выражение. Например, выражение

sin(2*12)/7+9^2

должно вычисляться по следующим шагам:
1.    Вычислить sin(2*12)/7
2.    Вычислить 9^2
3.    Затем применить оператор + к двум результатам
4.    Чтобы вычислить пункт 1, надо:
1.    Вычислить sin(2*12)
2.    Вычислить 7
3.    Затем поделить два результата
4.    Чтобы вычислить пункт 4.1, надо:
1.    Вычислить 2*12
2.    Применить функцию sin к результату
3.    Чтобы вычислить пункт 4.4.1, надо:
1.    Применить оператор * к двум операндам 2 и 12

Сканирующая функция принимает входное выражение и массив операндов в зависимости от их приоритета, как в следующем коде:

///////////////////////////////////////////////////////////
// GetOperator: сканирует входную строку,
// чтобы найти входные операторы
///////////////////////////////////////////////////////////
int GetOperator(IN LPCSTR lpcs, IN LPCSTR lpcsOperators[])
{
    for(int nIndex = 0; lpcsOperators[nIndex]; nIndex++)
    {
        int nOpen = 0;
        // сканируется выражение с конца
        LPSTR p = (LPSTR)lpcs+strlen(lpcs)-1;
        // цикл сообщает о достижении начала выражения
        while(p >= lpcs)
        {
            // проверка на закрывающую скобку
            if(*p == ')')
                nOpen++;
            // проверка на открывающую скобку
            else    if(*p == '(')
                nOpen--;
            // проверка на оператор
            else if(nOpen == 0 && strchr(lpcsOperators[nIndex], *p) != NULL)
                // проверка, является ли оператор меткой знака
                if((*p != '-' && *p != '+') ||
                   (p != lpcs && IsRightSign(*(p-1),
                         lpcsOperators, nIndex+1)))
                    // возвращается индекс оператора
                    return (int)(p-lpcs);
            p--;
        }
    }
    // оператор не найден
    return -1;
}

Очевидно, что это рекурсивная операция. Чтобы сделать это, надо сформулировать выражение в виде набора операторов и операндов и определенным образом обработать этот набор, чтобы получить конечный результат. Основные шаги для осуществления этого таковы:

•    Найти операторы
•    Разбить выражение по индексу оператора
•    Добавить два полученных операнда в набор
•    Если операторы не найдены, проверить, начинается ли выражение с имени функции

Применение предыдущих шагов дает следующий набор:

+

/

sin(2*12)

*

2

12

7

 

^

 

9

 

2

 

Для более подробного описания рассматривается код, разбирающий выражение и заполняющий набор:

bool FillStack(LPCSTR lpcsInput, vector<EXPRESSIONITEM*>& vStack)
{
    // массив операторов от верхнего до нижнего
    priority LPCSTR lpcsOperators[] = { "+-", "*/", "^%", NULL };

    // вставить первый входной элемент в набор
    vStack.push_back(new ExpressionItem(lpcsInput));
    // цикл проходит по набору выражения, чтобы проверить,
    // можно ли разбить выражение на два запроса
    for(int nIndex = 0; nIndex < (int)vStack.size(); nIndex++)
        // проверяет, является ли элемент выражения оператором
        if(vStack[nIndex]->m_cOperator == 0)
        {
            // копируется строка выражения
            CString str = vStack[nIndex]->m_strInput;
            // выражение разбирается, чтобы найти операторы
            int nOpIndex = GetOperator(str, lpcsOperators);
            if(nOpIndex != -1)
            {
                // выражение разбивается на
                // два запроса по индексу оператора
                vStack[nIndex]->m_cOperator = str[nOpIndex];
                // добавляется левый операнд оператора
                // как новое выражение
                vStack.insert(vStack.begin()+nIndex+1,
                  new ExpressionItem(str.Left(nOpIndex)));
                // добавляется правый операнд оператора
                // как новое выражение
                vStack.insert(vStack.begin()+nIndex+2,
                  new ExpressionItem(str.Mid(nOpIndex+1)));
            }
            else
            // проверяется, начинается ли строка выражения
            // с функции или круглых скобок
                if((vStack[nIndex]->m_nFunction = GetFunction(str,
                        vStack[nIndex]->m_nSign)) == 0
                        && vStack[nIndex]->m_nSign == 0)
                    // удаляются круглые скобки, и выражение повторно просматривается
                    vStack[nIndex--]->m_strInput =
                             str.Mid(1, str.GetLength()-2);
    }
    return true;
}

Это простейший код разбора математического выражения. Его комментарии собраны по пунктам ниже:

•    Вставить первый входной элемент в набор
•    Цикл проходит по набору выражения, чтобы проверить, можно ли разбить выражение на два запроса
•    Проверить, является ли элемент выражения оператором (+-*/^)
•    Разобрать выражение, чтобы найти операторы
•    Если оператор найден:
    o    Разбить выражение на два запроса по индексу оператора
    o    Добавить левый операнд оператора как новое выражение
    o    Добавить правый операнд оператора как новое выражение
•    Иначе
    o    Проверить, начинается ли строка выражения с функции или круглых скобок
    o    Если функция найдена, установить ее номер в данных выражения

Применение формул дифференцирования к набору

Применение формул дифференцирования –  значит выполнять итерации по набору, брать каждый оператор и применять правило дифференцирования (в зависимости от оператора)  к двум операндам оператора:

Выражение

Производная

u+v

du/dx+dv/dx

u-v

du/dx-dv/dx

u/v

(v*du/dx-u*dv/dx)/v^2

u*v

u*dv/dx+v*du/dx

c*u

c*du/dx

u^n

n*u^(n-1)*du/dx

c^u

c^u*ln(c)*du/dx

e^u

e^u*du/dx

Если любой операнд является оператором, то функция снова вызывается рекурсивно для дифференцирования всего выражения. Следующий код представляет эту функцию:

CString DifferentiateStack(vector<EXPRESSIONITEM*>& vStack, int& nExpression)
{
    ExpressionItem *pQI = vStack[nExpression++];
    if(pQI->m_cOperator)
    {
        // взять левый операнд
        CString u = vStack[nExpression]->GetInput();
        // взять производную от левого операнда
        CString du = DifferentiateStack(vStack, nExpression);
        // взять правый операнд
        CString v = vStack[nExpression]->GetInput();
        // взять производную от правого операнда
        CString dv = DifferentiateStack(vStack, nExpression);

        if(du == '0')    // u - константа
            ...
        else    if(dv == '0')    // v - константа
            ...
        else
            switch(pQI->m_cOperator)
            {
            case '-':    // d(u-v) = du-dv
            case '+':    // d(u+v) = du+dv
                pQI->m_strOutput = '('+du+pQI->m_cOperator+dv+')';
                break;
            case '*':    // d(u*v) = u*dv+du*v
                pQI->m_strOutput = '('+u+'*'+dv+'+'+du+'*'+v+')';
                break;
            case '/':    // d(u/v) = (du*v-u*dv)/v^2
                pQI->m_strOutput = '('+du+'*'+v+'-'+u+'*'+dv+")/("+v+")^2";
                break;
            case '^':    // d(u^v) = v*u^(v-1)*du+u^v*ln(u)*dv
                pQI->m_strOutput = '('+v+'*'+u+"^("+v+"-1)*"+ du+
                                      '+'+u+'^'+v+"*ln("+u+")*"+dv+')';
                break;
            }
    }
    else
        // взять производную выражения
        pQI->GetDifferentiation();
    // вернуть полученную производную
    return pQI->m_strOutput;
}

Оптимизация уравнения

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

•    Заменить "--" на "" или "+"
•    Заменить "+-" на "-"
•    Заменить "((....))" на "(....)"
•    Удалить все "1*"
•    Удалить все "*1"
•    Удалить все "^1"
•    Удалить все "1*"
•    Удалить ненужные круглые скобки

Псевдокод дифференцирования

Differentiate(input)
{
    Stack = FillStack(input)
    output = DifferentiateStack(Stack)
    Optimize(output)
    return output
}

FillStack(input)
{
    operators[] { "+-", "*/", "^%" }

    stack.push(input)
    loop( n = 1  to stack.size() )
    {
        if stack[n] is not operator
            if GetOperator(stack[n],
               operators) success
            {
                Split stack[n]
                stack.Insrt(left  operand)
                stack.Insrt(right operand)
            }
            else
                GetFunction(stack[n])
    }
}

DifferentiateStack(stack, index)
{
    if stack[index] is operator
    {
        index++
        u = stack[index]
        du = DifferentiateStack(stack, index)
        v = stack[index]
        dv = DifferentiateStack(stack, index)

        if operator = '-' or '+'
            output = du+operator+dv
        else    if operator = '*'
            output = u*dv+du*v
        else    if operator = '/'
            output = (du*v-u*dv)/v^2
        else    if operator = '^'
            output = v*u^(v-1)*du+u^v*ln(u)*dv
    }
    else
        output = stack[index++].GetDifferentiation()

    return output
}

void Optimize(str)
{
    replace "--" with "" or "+"
    replace "+-" with "-"
    replace "((....))"  with "(....)"
    remove any 1*
    remove any *1
    remove any exponent equal 1
    remove unneeded parentheses

    if str changed then
        Optimize(str)
}

ExpressionItem::GetDifferentiation()
{
    if Function then
    {
        arument = Argument(input);
        if function = SIN
            output = Differentiate(arument)*cos(arument)
        else    if function = COS
            output = Differentiate(arument)*(-sin(arument))
        else
        ...
        }
    }
    else
    {
        if input = "x"
            output = "1"
        else    if input = "-x"
            output = "-1"
        else    if input is numeric
            output = "0"
        else
            output = "d"+input+"/dx"
    }
}

Вкладки программы

Дифференцирование

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

e^sin(pi/3)/tan(x) оптимизируется до 2.377442/tan(x) перед дифференцированием, чтобы получить результат

(-2.377442*sec(x)^2)/(tan(x))^2.

Набор дифференцирования

Вычисление

Эта вкладка – символьный калькулятор, способный вычислить любое сложное математическое выражение. Вычисление содержит только две константы:

e = 2.7182818284590452353602874713527
pi = 3.1415926535897932384626433832795

Набор вычисления

Этот вид используется для просмотра набора деталей вычисления. c( ) означает "вычисление чего-либо". Содержимое набора расположено сверху вниз. Шаги вычисления являются шагами дифференцирования, за исключением формул операторов. Данный вид показывает части выражения, которые можно вычислить и которые нельзя вычислить, как в примере на рисунке.

Чертеж

Эта вкладка используется для просмотра кривой входной функции и кривой ее производной. Кривую входной функции можно просмотреть, нажав на "Чертеж". Вид содержит горизонтальные и вертикальные оси, а также линии сетки. Можно перейти на любую точку в виде, чтобы увидеть ее координаты в левом верхнем углу, и можно увеличить любую часть, выбрав ее мышью. Границы горизонтальных и вертикальных осей обновляются автоматически, чтобы включить точки входной функции в вид. Эти границы можно изменить, щелкнув правой кнопкой мыши по виду.

Справка

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


Интересные особенности

1.    Вычисление части выражения

Преимущество приведенного парсера выражения перед другими парсерами в интернете заключается в том, что он способен вычислить часть выражения и не трогать другую часть, включающую в себя любые переменные. Например, берется выражение sin(45+sin(2))/tan(x). Если ввести это выражение в любой бесплатный парсер в интернете, он даст ошибку, указывающую неправильный символ x. Однако приведенный парсер даст 0.937227/tan(x). Данный парсер пытается оптимизировать сложное выражение перед дифференцированием и вычисляет любую часть, поддающуюся вычислению. Итак, производная sin(45+sin(2))/tan(x) - (-0.937227*sec(x)^2)/(tan(x))^2.

2.    Приоритет оператора

Если в сложном выражении есть несколько операторов, приоритет оператора определяет последовательность, в которой выполняются операции. Порядок выполнения может сильно влиять на итоговое значение. Операторы имеют следующие уровни приоритета. Оператор более высокого уровня вычисляется раньше оператора более низкого уровня:

o    + (положительный), - (отрицательный)
o    ^ (показатель степени), % (остаток целочисленного деления)
o    * (умножение), / (деление)
o    + (сложение), - (вычитание)

Если два оператора в выражении имеют одинаковый уровень приоритета оператора, они вычисляются слева направо исходя из их положения в выражении. Например, в 4 - 2 + 27 оператор вычитания вычисляется раньше оператора сложения для получения 2 + 27, что дает результат выражения 29.

Круглые скобки применяются для отмены определенного приоритета операторов в выражении. Все внутри круглых скобок вычисляется первым, чтобы получить единственное значение, прежде чем это значение сможет использовать любой оператор вне круглых скобок. Например, 2 * 4 + 5 пересчитывается в 8 + 5, что дает результат выражения 13. Выражение 2 * (4 + 5) пересчитывается в 2 * 9; из-за скобок сложение выполняется первым, и результат выражения равен 18.

Если в выражении есть вложенные скобки, наиболее глубоко вложенное выражение вычисляется первым. Пример 2 * (4 + (5 - 3) ) содержит вложенные скобки, с выражением 5 – 3 в наиболее глубоко вложенной паре круглых скобок. Это выражение дает значение 2. Затем оператор сложения (+) прибавляет этот результат к 4, что дает значение 6. Наконец, 6 умножается на 2, давая результат выражения 12.

3.    Увеличение

В области чертежа можно выбрать любой участок для увеличения. Например, функция sin(100*x) имеет очень большую вариацию, как на следующем рисунке:

С помощью многократного увеличения можно увидеть более мелкие детали, как на рисунках ниже:

Замечание: Проверяйте разрешение горизонтальной оси после каждого увеличения.