Двумерная Паде-аппроксимация

[Fizykochemik]

Уважаемые коллеги!

Одномерная Паде-аппроксимация вещь подробно изученная и доступная во многих программах и библиотеках.
У меня возникла необходимость аппроксимировать двумерную функцию $f(x,y)$ в виде $\frac{P^{m_1,m_2}(x,y)}{Q^{n_1,n_2}(x,y)}$.
Для двумерного и более случаев я не нашел доступных программ, поэтому нужно писать самому.

Например, алгоритм для одномерного случая описан в Numerical Recipes (разд. 5.12), думаю обобщить его на двумерных случай.
Там главная соль в том что получающаяся СЛАУ чаще всего бывает сингулярная. Буду признателен за любые советы.
В моем случае $f(x,y)$ это полином по определению, если это имеет значение.

[Markiyan Hirnyk]

Результаты поиска в Гугле см. тут.
Последовательное применение одномерной аппроксимации Падэ дает удовлетворительные результаты. Это механизировано в Мэйпле, Математике и Матлабе. Вот пример, сделанный в Мэйпле:

with(numapprox):
A := pade(pade((sin(x)+cos(y))/(exp(x*y)+1), x, [2, 2]), y, [2, 2]);
A := (1/48)*(8064*x^10+77760*x^9+296352*x^8+751680*x^7+1518912*x^6+2363904*x^5+2861568*x^4+3732480*x^3+2985984*x^2+2239488*x+2239488+(-17856*x^10-71568*x^9-168480*x^8-415584*x^7-811296*x^6-964224*x^5-995328*x^4-715392*x^3-373248*x^2-186624*x)*y+(-18048*x^10-68616*x^9-73584*x^8-89424*x^7-203904*x^6-274752*x^5-554688*x^4-248832*x^3-1244160*x^2+186624*x-933120)*y^2)/(336*x^10+1224*x^9+5004*x^8+13392*x^7+27000*x^6+44064*x^5+72576*x^4+46656*x^3+124416*x^2+93312+(264*x^10+1122*x^9+2556*x^8+6876*x^7+6804*x^6+14256*x^5+2592*x^4+40176*x^3-7776*x^2+38880*x)*y+(208*x^10+951*x^9+1542*x^8+4284*x^7+10080*x^6+3456*x^5+29376*x^4+29808*x^2+7776)*y^2)

Не уверен, что результат является двухмерной аппроксимантой Падэ, но приближение хорошее

Optimization:-Maximize((A-(sin(x)+cos(y))/(exp(x*y)+1))^2, x = 0 .. (1/2)*Pi, y = 0 .. (1/2)*Pi);
 [0.0231303512590284588,    [x = HFloat(1.5707963267949), y = HFloat(1.5707963267949)]]

[Fizykochemik]

Спасибо за ответ! Думаю, как временное средство пригодиться, но для уверенности лучше сделать честную двумерную Паде аппр. и сравнить
с двухшаговым методом. В случае если $f(x,y)=f_1(x)+f_2(y)+g(x,y)$ для функций $f_1(x)$ и $f_2(y)$ можно использовать одномерную, а для $g(x,y)$ двумерную аппр.

Есть проблема, связанная с эмпиричностью выбора $m,n$, но так как я знаю поведение функции, то я могу выбрать подходящие $m,n$.
Насколько я понимаю, это проблема общая для Паде аппр.