Но интересно другое - то что нередко мои программы написанные лет 15 назад, остаются state-of-the-art - всё равно лучше по многим параметрам, чем почти всё что написано на эту тему на Матлабе с тех пор. На днях получил е-мейл насчёт программы для численного интегрирования, которую я запостил в матлабовскую ньюсгруппу ещё в 1995, потом кто-то в 1999 залил на сайт Mathworks:
Dear Dr. Kirill Pankratov,
I am a PhD student at MIT .
I was looking at your quads.m and quadr.m files online from
http://www.mathworks.com/matlabcentral/linkexchange/links/41-quadrmquadsm
I get the following error message when I try to run then, I am not sure why they are not working.
[Q,xo,fo,s] = quads('1/sqrt(x)/(x-2)/(x-4)^2',0,5);
??? Subscript indices must either be real positive integers or logicals.
Error in ==> lininsrt at 75
xo(mask,:) = x;
Error in ==> quads at 160
[x_o,mask] = lininsrt(x_o,F);
Is it because you wrote these functions in 1999 and now I am using a newer Matlab version?
I really appreciate your help,
Thank you very much,
G*******
Чуть-чуть покопался, изменил что надо - работает. Послал просившему, тот оказался вполне удовлетворён.
Заодно посмотрел как она ведёт себя по сравнению с тем что в Матлабе имеется сейчас. Выяснилось, что, хотя за эти годы программы численного интегрирования улучшились, они недотягивают до того что я написал 15 лет назад.
Мои программы были особенно полезны для "трудных" функций - в том числе с высокочастотными колебаниями, очень широким динамическим диапазоном и сингулярностями.
Как правило мои программа одновременно быстрее, точнее и надёжнее чем собственная матлабовская
Например, моя функция quadr запросто может проинтегрировать 10 тысяч синусоид (точный ответ - 2.0):
>> quadr('sin(x)',0,10001*pi,1e-6)
Recursion level limit reached. Singularity likely.
ans =
2.0001
В то время как родную матлабовскую quad заносит на виражах в бесконечность даже на 1/25 этого диапазона:
>> quad('sin(x)',0,401*pi,1e-6)
Warning: Maximum function count exceeded; singularity likely.
> In quad at 92
ans =
-105.5199
Моя функция может даже автоматически найти, проанализировать и, если возможно, проинтегрировать сингулярности. Например, она может определить является ли сингулярность интегрируемой, интегрируемой в смысле главного значения (типа 1/х) или неинтегрируемой. Ничего подобного матлабовские родные программы не делают даже близко.
Из описания прилагаемого к моим программам:
>> [Q,xo,fo,s] = quads('1/x',-1,2);
Singularity at x = -0.000000, integrable in a principal value sense
Check the accuracy:
>> Q/log(2)-1
ans =
-1.1324e-06
(LOG(2) is a value of integral from 1 to 2. Integral from -1
to 1 is zero in a principal value sense).
...
Let's try a function with several singularities of different
types:
>> [Q,xo,fo,s] = quads('1/sqrt(x)/(x-2)/(x-4)^2',0,5);
Singularity at x = 0.000008, integrable
Singularity at x = 2.000000, integrable in a principal value sense
Singularity at x = 4.000000, non-integrable
Now try a function which is not strictly polynomial or
rational:
>> [Q,xo,fo,s] = quads('x*tan(x)',0,pi*4);
Singularity at x = 1.570796, integrable in a principal value sense
Singularity at x = 4.712389, integrable in a principal value sense
Singularity at x = 7.853982, integrable in a principal value sense
Singularity at x = 10.995574, integrable in a principal value sense
One can easily check that singularities are at the points
pi/2, 3*pi/2, 5*pi/2, 7*pi/2:
>> s(1,:)/pi
ans =
0.5000 1.5000 2.5000 3.5000
По-моему, круто! :)