Pull to refresh
28
0
Петр Калинин @pkalinin

Физик-программист и учитель

Send message
Вы, похоже, основным достоинством символьных вычисления полагаете, что для них есть готовые пакеты. Для автоматического дифференцирования готовые пакеты тоже есть, я про это неоднократно говорю в тексте и указываю, где их можно взять. И если вы думаете, что я предлагаю всегда использовать «написанный на коленке код Derivable», то вы ошибаетесь. Ясно, что в определенных приложениях лучше взять готовый и более надежный код — но вот вам пожалуйста, пройдите по паре ссылок из поста и найдете их, готовых, огромное количество.

Но мне кажется столь же ясным, что во многих приложениях проще написать код на 50 строчек и оттестировать его, чем подключать стороннюю библиотеку. На тестирование того, что вы правильно понимаете стороннюю библиотеку и на тестирование кода связки у вас уйдет не меньше времени, чем на тестирование корректности кода Derivable. И если ваша задача — считать производную для функции, записанной несколькими операторами присваивания с, возможно, использованием вспомогательных функций (типичный пример в комменте), то это будет работать, и не потребует (что, на мой взгляд, весьма важно) никакого изменения в промежуточный код, кроме замены double на Derivable.
Поправка: в коде таки баг, надо как минимум
    xd = Derivable::IndependendVariable(x);
    Derivable u = foo(xd, Q * xd / 2);
    Derivable v = bar(xd / t, constants::c); // скорость света

А вообще, mayorovp выше написал еще правильнее.
Более того (хотя Mrrl привел тоже хороший пример), пусть вы, например, решаете какое-то уравнение f(x,a)=0 относительно x тем же методом Ньютона. И хотите получать не только корень, но и его производную по a. Символьно вы это никак не сделаете, т.к. условие выхода из цикла зависит от текущего (на итерации) значения x. (Ну или будете вынуждены фиксировать количество итераций). А с Derivable все просто:
// решает уравнение f(x,a)=x*x-a=0 и возвращает корень и его производную по a
Derivable sqrtWithDeriv(double a) { 
    static const double EPS = 1e-6;
    Derivable ad = Derivable::IndependendVariable(a);
    Derivable x = 1;
    while (1) {
        Derivable thisF = x*x - ad;
        if (fabs(thisF.getValue()) < EPS) return x;
        x = (ad + x * x)/2/x;
    }
}

(В этом примере производную df/dx я посчитал руками, но можно было и сделать Derivable по двум переменным и написать цикл так же, как и в примере.)
Реже придется делать перегруженную версию (как, например, для sin и cos).

Для элементарных функций придется для всех. Если же у вас есть доступ к исходному коду функции, то можно внутри нее все поменять.
Если вы про раздел 2.3.2, то там, насколько я понял, речь про символьное дифференцирование — запишем выражение строкой (ну или списком), и будем работать с этой строкой.

Автоматическое дифференцирование — несколько другой подход, и, пожалуй, ответ на ваш второй вопрос довольно хорошо иллюстрирует эту разницу. Как решить исходную задачу? Как и написано в посте,
главное везде все промежуточные переменные (и в том числе результаты промежуточных функций) сделать Derivable, остальной код даже менять не потребуется!

Пример (очень условный и выдуманный мною только что; компилируемость и работоспособность не проверял, но идея, думаю, понятна): пусть был код
double foo(double x, double y) {
    if (x < y) return x * x - y;
    else return y * y - x;
}

double bar(double x, double c) {
    return sin(x / c) * c / x;
}

double f(double x, double Q, double t) {
    double u = foo(x, Q * x / 2);
    double v = bar(x / t, constant::c); // скорость света
    return u / v;
}

и нам надо получить df/dx.

Заменяем везде, где речь идет про вычисляемые функции, double на Derivable:
Derivable foo(Derivable x, Derivable y) {
    if (x < y) return x * x - y;
    else return y * y - x;
}

Derivable bar(Derivable x, double c) { // считаем, что тут c никогда не зависит от x
    return sin(x / c) * c / x;
}

Derivable f(double x, double Q, double t) {
    Derivable u = foo(x, Q * x / 2);
    Derivable v = bar(x / t, constants::c); // скорость света
    return u / v;
}

Ну еще придется или перегрузить operator<, или сделать приведение типов от Derivable обратно к double.

Information

Rating
Does not participate
Location
Нижний Новгород, Нижегородская обл., Россия
Date of birth
Registered
Activity