Зависимость времени инт-ия от формы коэф. метода Рунге-Кутта

[Artem_]

Странная зависимость наблюдалась мной при численном интегрировании никоторой системы уравнений методом Рунге-Кутта. Если взять коэффициенты схемы в форме рациональных дробей, то общее время оказывалось Х с. Если те же самые коэффициенты в той же самой схеме подставить в вещественной форме (с плавающей точкой), то время интегрирования при всех прочих равных условий заняло 90*Х с.
Использовался метод Рунге-Кутта 4-5 порядка с адаптацией шага.
Я думаю, что и в методах более низкого порядка будет наблюдаться подобный разрыв.

А ведь это странно! Ведь в какой бы форме мы не подставляли бы коэффициенты при расчете они все равно интерпретатором (системы программирования или системы компьютерной математики) в конечном итоге преобразуются в форму с плавающей запятой, разве нет?

[Artem_]

Попросили указать среду. Среда Matlab 7.0.1.24704 (R14) (SP1)

[photon]

Ага, тогда давайте код - посмотрим
Вы использовали функцию ode45()?

[Artem_]

Нет, конечно, я использовал не функцию ode45, поскольку речь идет о КОЭФФИЦИЕНТАХ алгоритма, которые черт знает где в стандартной функции ode45 прописаны. Я писал стандартные алгоритмы 23, 34, 45 с нуля в m-файле.

[photon]

1) Это я догадался. Я имел в виду: Вы не пробовали использовать стандартную и сравнить с ней свои? Скорее всего, оба Ваших варианта проиграют.

Вопрос: Зачем Вам писать свои, если есть готовое? (Конечно, если надо, то надо)

2) Давайте коды m-файлов: посмотрим, поэкспериментируем вместе.

А вот Вам цитата из файла ode45.m, начиная со строки 244:

% Initialize method parameters.
pow = 1/5;
A = [1/5; 3/10; 4/5; 8/9; 1; 1];
B = [
    1/5         3/40    44/45   19372/6561      9017/3168       35/384
    0           9/40    -56/15  -25360/2187     -355/33         0
    0           0       32/9    64448/6561      46732/5247      500/1113
    0           0       0       -212/729        49/176          125/192
    0           0       0       0               -5103/18656     -2187/6784
    0           0       0       0               0               11/84
    0           0       0       0               0               0
    ];
E = [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40];

[незваный гость]

Artem_ писал(а):

Ведь в какой бы форме мы не подставляли бы коэффициенты при расчете они все равно интерпретатором (системы программирования или системы компьютерной математики) в конечном итоге преобразуются в форму с плавающей запятой, разве нет?

Если говорить о представления данных в компьютере, то разница существует, причем на нескольких уровнях. $1/5$ -- точное число, в то время как $0.2$ уже будет представлено приближенно. Кроме того, на уровне команд / интерпертатора это могут быть очень разные операции. Ошибки же округления могут вести к уменьшению шага интегрирования, что, в свою очередь, к замедлению.

Ну и наконец, самый простой вариант -- проверить, все ли коэффициенты подсчитаны правильно. Ошибка в вычислении коэффициентов (или их представление с недостаточной точностью) -- самое простое объяснение. Был у меня случай, один товарищ $\pi$ тройкой аппроксимировал... Долгонько искали...

[photon]

незванный гость писал(а):

:evil:

Artem_ писал(а):

Ведь в какой бы форме мы не подставляли бы коэффициенты при расчете они все равно интерпретатором (системы программирования или системы компьютерной математики) в конечном итоге преобразуются в форму с плавающей запятой, разве нет?

Если говорить о представления данных в компьютере, то разница существует, причем на нескольких уровнях. $1/5$ -- точное число, в то время как $0.2$ уже будет представлено приближенно. Кроме того, на уровне команд / интерпертатора это могут быть очень разные операции. Ошибки же округления могут вести к уменьшению шага интегрирования, что, в свою очередь, к замедлению.

На самом деле, ситуация действительно странная: у MatLAB по умолчанию все числа имеют тип double, т.е., по идее, дополнительная операция деления должна замедлять работу. Я сейчас попробовал на своем компьютере два простеньких примера:

tic, for i=1:1e5 a=i*0.2;end, toc

и

tic, for i=1:1e5 a=i*(1/5);end, toc

В первом случае устойчиво получается быстрее примерно в 1.3 раза, а вот Artem_ говорит, что у него с дробями быстрее, причем гораздо, и в стандартном MatLAB-файле дроби, не думаю, что тут разработчики прошляпили. В общем, действительно странно.

[Comga]

Artem_ писал(а):

Странная зависимость наблюдалась мной при численном интегрировании никоторой системы уравнений методом Рунге-Кутта. Если взять коэффициенты схемы в форме рациональных дробей, то общее время оказывалось Х с. Если те же самые коэффициенты в той же самой схеме подставить в вещественной форме (с плавающей точкой), то время интегрирования при всех прочих равных условий заняло 90*Х с.
Использовался метод Рунге-Кутта 4-5 порядка с адаптацией шага.

1) А результаты-то совпадают?
2) Раз уж все пишете сами, перепишите на Фортран (С), должно работать еще быстрее
3) Как известно, программ (речь в данном случае о matlab) без ошибок не бывает
4) Привели бы оба текста - проще было бы советовать

[Artem_]

Огромное спасибо г-ну photon_у за ценный совет.

Действительно, чтобы не грешить на ошибки, допущенные в своем коде, давайте изменим код стандартной Matlab_овской функции и посмотрим что получиться.
Для этого:
1 копируем из каталога
\MATLAB701\toolbox\matlab\funfun\ в свой рабочий каталог файл ode45.m переименовав его, например, в ode451.m.
2 исправляем в этом файле, начиная со строки 265, 4 вектора: pow, A, B, E заменив в них дроби на числа с плавающей точкой.
3 копируем из каталога
\MATLAB701\toolbox\matlab\funfun\private\ в свой рабочий каталог функции odearguments.m, ntrp45.m, odemass.m, odeevents.m, odefinalize.m, поскольку ode45 обращается к этим функциям в процессе.
Затем вы можете попробовать решить нек. задачу решателем ode45 и с теми же условиями решателем ode451.
Даже не используя циклы вы увидите разницу. Во-первых, итоговый вектор решения будет длиннее при использовании ode451, во-вторых, время расчета увеличиться.
Увеличение времени при одном расчете вашей задачи будет возможно несущественным. В моем случае оно составляло от 10% до 20%, но если Вам нужно (как в моем случае) решать систему в цикле или на большом временном интервале, то получается ого-го.

Так вот и получается, что наша небольшая невинная шутка - замена дробей на плавующую точку приводит в Matlab_е к неприятностям.

Ваши комментарии?