Using Maxima to solve certain mathematical problems. (Russian)

Table 1: Clock summary at [2023-08-27 Sun 23:42]
Headline	Time
Total time	13:04
Решение	13:04
Найти 20 ортогональных полиномов…		11:25
Используя формулу Родрига			3:55
Используя рекуррентную формулу			0:47
Используя производящую функцию			1:27
Используя гипергеометрическую функцию…			5:16
Сделать трансляцию на язык Си…		1:39

В этом файле я бы хотел решить несколько задач по математике с использованием системы символьных вычислений Maxima. Задачи взяты из учебника "Ильина-Силаев".

2. Решение

2.1. Найти 20 ортогональных полиномов (0-19) Лежандра

2.1.1. Используя формулу Родрига

Формула Родрига: https://en.wikipedia.org/wiki/Rodrigues%27_formula

\[ \frac{(-1)^{n}}{2^{n}\cdot n!} \frac{\partial}{dx}(1 - x^{2})^{n} \]

P_rodrigues[n]:= ev(((-1)^n)/((2^n) * (n!)) * diff((1 - x^2)^n, x, n), diff, expand);

P_rodrigues[n]:= ev(((-1)^n)/((2^n) * (n!)) * diff((1 - x^2)^n, x, n), diff, expand);
for i:0 thru 19 do ( simp:true, fr:rat(P_rodrigues[i]),
		     simp:false, tex((1/denom(fr)) * num(fr)));

\[ {{1}\over{1}}\,1 \] \[ {{1}\over{1}}\,x \] \[ {{1}\over{2}}\,\left(3\,x^2-1\right) \] \[ {{1}\over{2}}\,\left(5\,x^3-3\,x\right) \] \[ {{1}\over{8}}\,\left(35\,x^4-30\,x^2+3\right) \] \[ {{1}\over{8}}\,\left(63\,x^5-70\,x^3+15\,x\right) \] \[ {{1}\over{16}}\,\left(231\,x^6-315\,x^4+105\,x^2-5\right) \] \[ {{1}\over{16}}\,\left(429\,x^7-693\,x^5+315\,x^3-35\,x\right) \] \[ {{1}\over{128}}\,\left(6435\,x^8-12012\,x^6+6930\,x^4-1260\,x^2+35 \right) \] \[ {{1}\over{128}}\,\left(12155\,x^9-25740\,x^7+18018\,x^5-4620\,x^3+ 315\,x\right) \] \[ {{1}\over{256}}\,\left(46189\,x^{10}-109395\,x^8+90090\,x^6-30030 \,x^4+3465\,x^2-63\right) \] \[ {{1}\over{256}}\,\left(88179\,x^{11}-230945\,x^9+218790\,x^7-90090 \,x^5+15015\,x^3-693\,x\right) \] \[ {{1}\over{1024}}\,\left(676039\,x^{12}-1939938\,x^{10}+2078505\,x^ 8-1021020\,x^6+225225\,x^4-18018\,x^2+231\right) \] \[ {{1}\over{1024}}\,\left(1300075\,x^{13}-4056234\,x^{11}+4849845\,x ^9-2771340\,x^7+765765\,x^5-90090\,x^3+3003\,x\right) \] \[ {{1}\over{2048}}\,\left(5014575\,x^{14}-16900975\,x^{12}+22309287 \,x^{10}-14549535\,x^8+4849845\,x^6-765765\,x^4+45045\,x^2-429 \right) \] \[ {{1}\over{2048}}\,\left(9694845\,x^{15}-35102025\,x^{13}+50702925 \,x^{11}-37182145\,x^9+14549535\,x^7-2909907\,x^5+255255\,x^3-6435\, x\right) \] Unable to evaluate predicate 1 + (1 + (1

(1 + (1 + (1 + (1 + (1 + (1 + (1 + (1
(1 + (1 + (1 + (1 + (1 + 0))))))))))))))) > 19

– an error. To debug this try: debugmode(true);

Очень любопытно, что на степени 16 случилось переполнение непонятно чего.

P_rodrigues[n]:= ev(((-1)^n)/((2^n) * (n!)) * diff((1 - x^2)^n, x, n), diff, expand);
  for i:0 thru 19 do tex(P_rodrigues[i]);

\[ 1 \] \[ x \] \[ {{3\,x^2}\over{2}}-{{1}\over{2}} \] \[ {{5\,x^3}\over{2}}-{{3\,x}\over{2}} \] \[ {{35\,x^4}\over{8}}-{{15\,x^2}\over{4}}+{{3}\over{8}} \] \[ {{63\,x^5}\over{8}}-{{35\,x^3}\over{4}}+{{15\,x}\over{8}} \] \[ {{231\,x^6}\over{16}}-{{315\,x^4}\over{16}}+{{105\,x^2}\over{16}}- {{5}\over{16}} \] \[ {{429\,x^7}\over{16}}-{{693\,x^5}\over{16}}+{{315\,x^3}\over{16}}- {{35\,x}\over{16}} \] \[ {{6435\,x^8}\over{128}}-{{3003\,x^6}\over{32}}+{{3465\,x^4}\over{ 64}}-{{315\,x^2}\over{32}}+{{35}\over{128}} \] \[ {{12155\,x^9}\over{128}}-{{6435\,x^7}\over{32}}+{{9009\,x^5}\over{ 64}}-{{1155\,x^3}\over{32}}+{{315\,x}\over{128}} \] \[ {{46189\,x^{10}}\over{256}}-{{109395\,x^8}\over{256}}+{{45045\,x^6 }\over{128}}-{{15015\,x^4}\over{128}}+{{3465\,x^2}\over{256}}-{{63 }\over{256}} \] \[ {{88179\,x^{11}}\over{256}}-{{230945\,x^9}\over{256}}+{{109395\,x^ 7}\over{128}}-{{45045\,x^5}\over{128}}+{{15015\,x^3}\over{256}}-{{ 693\,x}\over{256}} \] \[ {{676039\,x^{12}}\over{1024}}-{{969969\,x^{10}}\over{512}}+{{ 2078505\,x^8}\over{1024}}-{{255255\,x^6}\over{256}}+{{225225\,x^4 }\over{1024}}-{{9009\,x^2}\over{512}}+{{231}\over{1024}} \] \[ {{1300075\,x^{13}}\over{1024}}-{{2028117\,x^{11}}\over{512}}+{{ 4849845\,x^9}\over{1024}}-{{692835\,x^7}\over{256}}+{{765765\,x^5 }\over{1024}}-{{45045\,x^3}\over{512}}+{{3003\,x}\over{1024}} \] \[ {{5014575\,x^{14}}\over{2048}}-{{16900975\,x^{12}}\over{2048}}+{{ 22309287\,x^{10}}\over{2048}}-{{14549535\,x^8}\over{2048}}+{{4849845 \,x^6}\over{2048}}-{{765765\,x^4}\over{2048}}+{{45045\,x^2}\over{ 2048}}-{{429}\over{2048}} \] \[ {{9694845\,x^{15}}\over{2048}}-{{35102025\,x^{13}}\over{2048}}+{{ 50702925\,x^{11}}\over{2048}}-{{37182145\,x^9}\over{2048}}+{{ 14549535\,x^7}\over{2048}}-{{2909907\,x^5}\over{2048}}+{{255255\,x^3 }\over{2048}}-{{6435\,x}\over{2048}} \] \[ {{300540195\,x^{16}}\over{32768}}-{{145422675\,x^{14}}\over{4096}} +{{456326325\,x^{12}}\over{8192}}-{{185910725\,x^{10}}\over{4096}}+ {{334639305\,x^8}\over{16384}}-{{20369349\,x^6}\over{4096}}+{{ 4849845\,x^4}\over{8192}}-{{109395\,x^2}\over{4096}}+{{6435}\over{ 32768}} \] \[ {{583401555\,x^{17}}\over{32768}}-{{300540195\,x^{15}}\over{4096}} +{{1017958725\,x^{13}}\over{8192}}-{{456326325\,x^{11}}\over{4096}}+ {{929553625\,x^9}\over{16384}}-{{66927861\,x^7}\over{4096}}+{{ 20369349\,x^5}\over{8192}}-{{692835\,x^3}\over{4096}}+{{109395\,x }\over{32768}} \] \[ {{2268783825\,x^{18}}\over{65536}}-{{9917826435\,x^{16}}\over{ 65536}}+{{4508102925\,x^{14}}\over{16384}}-{{4411154475\,x^{12} }\over{16384}}+{{5019589575\,x^{10}}\over{32768}}-{{1673196525\,x^8 }\over{32768}}+{{156165009\,x^6}\over{16384}}-{{14549535\,x^4}\over{ 16384}}+{{2078505\,x^2}\over{65536}}-{{12155}\over{65536}} \] \[ {{4418157975\,x^{19}}\over{65536}}-{{20419054425\,x^{17}}\over{ 65536}}+{{9917826435\,x^{15}}\over{16384}}-{{10518906825\,x^{13} }\over{16384}}+{{13233463425\,x^{11}}\over{32768}}-{{5019589575\,x^9 }\over{32768}}+{{557732175\,x^7}\over{16384}}-{{66927861\,x^5}\over{ 16384}}+{{14549535\,x^3}\over{65536}}-{{230945\,x}\over{65536}} \]

2.1.2. Используя рекуррентную формулу

Рекуррентная формула, взятая из того же самого учебного пособия Ильина-Силаев:

\[(n+1)P_{n+1} = (2n + 1)xP_{n}-nP_{n-1}\]

Рекурсия неплохо работает с мемоизацией, по-идее. Давайте попробуем помемоизировать.

\[n P_{n} = (2n-1)xP_{n-1}-(n-1)P_{n-2}\]

P_recurrent[n] := expand(((2*n - 1)*x*P_recurrent[n-1] - (n-1)*P_recurrent[n-2])/(n));
P_recurrent[0] : 1;
P_recurrent[1] : x;

P_recurrent[n] := expand(((2*n - 1)*x*P_recurrent[n-1] - (n-1)*P_recurrent[n-2])/(n));
P_recurrent[0] : 1;
P_recurrent[1] : x;
for i:0 thru 19 do tex(P_recurrent[i]);

\[ 1 \] \[ x \] \[ {{3\,x^2-1}\over{2}} \] \[ {{5\,x^3-3\,x}\over{2}} \] \[ {{35\,x^4-30\,x^2+3}\over{8}} \] \[ {{63\,x^5-70\,x^3+15\,x}\over{8}} \] \[ {{231\,x^6-315\,x^4+105\,x^2-5}\over{16}} \] \[ {{429\,x^7-693\,x^5+315\,x^3-35\,x}\over{16}} \] \[ {{6435\,x^8-12012\,x^6+6930\,x^4-1260\,x^2+35}\over{128}} \] \[ {{12155\,x^9-25740\,x^7+18018\,x^5-4620\,x^3+315\,x}\over{128}} \] \[ {{46189\,x^{10}-109395\,x^8+90090\,x^6-30030\,x^4+3465\,x^2-63 }\over{256}} \] \[ {{88179\,x^{11}-230945\,x^9+218790\,x^7-90090\,x^5+15015\,x^3-693 \,x}\over{256}} \] \[ {{676039\,x^{12}-1939938\,x^{10}+2078505\,x^8-1021020\,x^6+225225 \,x^4-18018\,x^2+231}\over{1024}} \] \[ {{1300075\,x^{13}-4056234\,x^{11}+4849845\,x^9-2771340\,x^7+765765 \,x^5-90090\,x^3+3003\,x}\over{1024}} \] \[ {{5014575\,x^{14}-16900975\,x^{12}+22309287\,x^{10}-14549535\,x^8+ 4849845\,x^6-765765\,x^4+45045\,x^2-429}\over{2048}} \] \[ {{9694845\,x^{15}-35102025\,x^{13}+50702925\,x^{11}-37182145\,x^9+ 14549535\,x^7-2909907\,x^5+255255\,x^3-6435\,x}\over{2048}} \] \[ {{300540195\,x^{16}-1163381400\,x^{14}+1825305300\,x^{12}- 1487285800\,x^{10}+669278610\,x^8-162954792\,x^6+19399380\,x^4- 875160\,x^2+6435}\over{32768}} \] \[ {{583401555\,x^{17}-2404321560\,x^{15}+4071834900\,x^{13}- 3650610600\,x^{11}+1859107250\,x^9-535422888\,x^7+81477396\,x^5- 5542680\,x^3+109395\,x}\over{32768}} \] \[ {{2268783825\,x^{18}-9917826435\,x^{16}+18032411700\,x^{14}- 17644617900\,x^{12}+10039179150\,x^{10}-3346393050\,x^8+624660036\,x ^6-58198140\,x^4+2078505\,x^2-12155}\over{65536}} \] \[ {{4418157975\,x^{19}-20419054425\,x^{17}+39671305740\,x^{15}- 42075627300\,x^{13}+26466926850\,x^{11}-10039179150\,x^9+2230928700 \,x^7-267711444\,x^5+14549535\,x^3-230945\,x}\over{65536}} \]

И всё?

2.1.3. Используя производящую функцию

Производящая функция – это вот эта функция: https://en.wikipedia.org/wiki/Generating_function

На википедии есть прямо инструкция о том, как её использовать для того, чтобы вычислить полиномы Лежандра:

https://en.wikipedia.org/wiki/Legendre_polynomials#Definition_via_generating_function

\[ \frac{1}{\sqrt{1-2xt+t^2}} = \sum_{n=0}^\infty P_n(x) t^n \]

P_generating[n]:= expand(part(taylor((1)/(sqrt(1 - 2*x*t + t^2)),t,0,19), [i+1])/(t^(i+1-1))) ;

P_generating[n]:= expand(part(taylor((1)/(sqrt(1 - 2*x*t + t^2)),t,0,19), [i+1])/(t^(i+1-1))) ;
for i:0 thru 19 do tex(P_generating[i]);

2.1.4. Используя гипергеометрическую функцию или вырожденную гипергеометрическую функцию

Гипергеометрическая функция – это очень классная функция, см https://en.wikipedia.org/wiki/Hypergeometric_function и https://dlmf.nist.gov/15 .

\[ F(a,b,c,z) = \sum_{n=0}^{\infty} \frac{(a)_{n}(b)_{n}}{(c)_{n}} \frac{z^{n}}{n!} \]

где $(a)_{n}$ – символ Похгаммера

\[(a)_{n}=\prod_{i=0}^{n-1}(a+i) \]

Выражение через неё полиномов Лежандра делается так:

\[ P_{n}= F\left(-n,n+1,1,\frac{1-x}{2}\right) \]

Радует, что в Maxima есть функция pochhammer.

Попробуем для начала реализовать гипергеометрическую функцию с нуля?

showtime: true;
load("simplify_sum");
F: sum( ((pochhammer(a,n)*pochhammer(b,n))/(pochhammer(c,n)))*((z^n)/(n!)) , n , 0 , inf );
tex(F);
P: subst([a=-p,b=p+1,c=1,z=((1-x)/2)], F);
simplify_sum(subst([p=1],P));

\[ \sum_{n=0}^{\infty }{{{\left(a\right)_{n}\,\left(b\right)_{n}\,z^{ n}}\over{\left(c\right)_{n}\,n!}}} \] (- 1) (2) (x - 1) n + 1 n + 1

----------------------–— non-rational term ratio to nusum 2 2 (- 1) (2) (n + 1) n n

Evaluation took 0.9851 seconds (1.0134 elapsed) using 36.178 MB.

Что-то результат так себе. Попробуем взять готовую гипергеометрическую функцию.

P_hypergeometric[n] :=  ev(hypergeometric([-n,n+1],[1],((1-x)/2)), diff, expand, simp);

P_hypergeometric[n] :=  ev(hypergeometric([-n,n+1],[1],((1-x)/2)), diff, expand, simp);
for i:0 thru 19 do tex(P_hypergeometric[i]);

Второй подход оказался намного быстрее.

2.1.5. Сравнить ответы

Okay, basically I will solve it by:

Partition #src blocks between implementation and tex print.
find some maxima way to compare different expressions.
make a for loop with expanding noweb links to polynomial definitions

I1: elapsed_real_time();
P_rodrigues[n]:= ev(((-1)^n)/((2^n) * (n!)) * diff((1 - x^2)^n, x, n), diff, expand);
I2: elapsed_real_time();
P_recurrent[n] := expand(((2*n - 1)*x*P_recurrent[n-1] - (n-1)*P_recurrent[n-2])/(n));
P_recurrent[0] : 1;
P_recurrent[1] : x;
I3: elapsed_real_time();
P_generating[n]:= expand(part(taylor((1)/(sqrt(1 - 2*x*t + t^2)),t,0,19), [i+1])/(t^(i+1-1))) ;
I4: elapsed_real_time();
P_hypergeometric[n] :=  ev(hypergeometric([-n,n+1],[1],((1-x)/2)), diff, expand, simp);
I5: elapsed_real_time();
for i:0 thru 19 do block([],
  /* tex(P_rodrigues[i]),*/
  /* tex(P_recurrent[i]), */
  ldisplay( rat(P_recurrent[i] - P_rodrigues[i]) = 0 ),
  /* tex(P_generating[i]), */
  ldisplay( rat(P_generating[i] - P_recurrent[i]) = 0 ),
  /* tex(P_hypergeometric[i]), */
  ldisplay( rat(P_hypergeometric[i] - P_generating[i]) = 0 ),
  printf(true, "~%----------------------~%"));

(%t1)/R/ 0 = 0 (%t2)/R/ 0 = 0 (%t3)/R/ 0 = 0

(%t4)/R/ 0 = 0 (%t5)/R/ 0 = 0 (%t6)/R/ 0 = 0

(%t7)/R/ 0 = 0 (%t8)/R/ 0 = 0 (%t9)/R/ 0 = 0

(%t10)/R/ 0 = 0 (%t11)/R/ 0 = 0 (%t12)/R/ 0 = 0

(%t13)/R/ 0 = 0 (%t14)/R/ 0 = 0 (%t15)/R/ 0 = 0

(%t16)/R/ 0 = 0 (%t17)/R/ 0 = 0 (%t18)/R/ 0 = 0

(%t19)/R/ 0 = 0 (%t20)/R/ 0 = 0 (%t21)/R/ 0 = 0

(%t22)/R/ 0 = 0 (%t23)/R/ 0 = 0 (%t24)/R/ 0 = 0

(%t25)/R/ 0 = 0 (%t26)/R/ 0 = 0 (%t27)/R/ 0 = 0

(%t28)/R/ 0 = 0 (%t29)/R/ 0 = 0 (%t30)/R/ 0 = 0

(%t31)/R/ 0 = 0 (%t32)/R/ 0 = 0 (%t33)/R/ 0 = 0

(%t34)/R/ 0 = 0 (%t35)/R/ 0 = 0 (%t36)/R/ 0 = 0

(%t37)/R/ 0 = 0 (%t38)/R/ 0 = 0 (%t39)/R/ 0 = 0

(%t40)/R/ 0 = 0 (%t41)/R/ 0 = 0 (%t42)/R/ 0 = 0

(%t43)/R/ 0 = 0 (%t44)/R/ 0 = 0 (%t45)/R/ 0 = 0

(%t46)/R/ 0 = 0 (%t47)/R/ 0 = 0 (%t48)/R/ 0 = 0

(%t49)/R/ 0 = 0 (%t50)/R/ 0 = 0 (%t51)/R/ 0 = 0

(%t52)/R/ 0 = 0 (%t53)/R/ 0 = 0 (%t54)/R/ 0 = 0

(%t55)/R/ 0 = 0 (%t56)/R/ 0 = 0 (%t57)/R/ 0 = 0

(%t58)/R/ 0 = 0 (%t59)/R/ 0 = 0 (%t60)/R/ 0 = 0

Кажется, всё сходится.

2.1.6. Собрать скорость работы в таблицу

Do the same loop as above
Measure runtime with elapsed_real_time()
* Maybe translate or compile maxima functions for speed.

  I1: elapsed_real_time();
  P_rodrigues[n]:= ev(((-1)^n)/((2^n) * (n!)) * diff((1 - x^2)^n, x, n), diff, expand);
  I2: elapsed_real_time();
  P_recurrent[n] := expand(((2*n - 1)*x*P_recurrent[n-1] - (n-1)*P_recurrent[n-2])/(n));
  P_recurrent[0] : 1;
  P_recurrent[1] : x;
  I3: elapsed_real_time();
  P_generating[n]:= expand(part(taylor((1)/(sqrt(1 - 2*x*t + t^2)),t,0,19), [i+1])/(t^(i+1-1))) ;
  I4: elapsed_real_time();
  P_hypergeometric[n] :=  ev(hypergeometric([-n,n+1],[1],((1-x)/2)), diff, expand, simp);
  I5: elapsed_real_time();
  printf(true, "rodrigues__time=~f~%", float(I2-I1));
  printf(true, "recurrent__time=~f~%", float(I3)-float(I2));
  printf(true, "generating_time=~f~%", float(I4)-float(I3));
  printf(true, "hypergeome_time=~f~%", float(I5)-float(I4));
  measure_time(expression) ::= block( [time: elapsed_real_time(), dummy:0],
    dummy:ev(expression) ,
    elapsed_real_time() - time );
  results: matrix( ['rodrigues], ['recurrent], ['generating], ['hypergeometric]);
  fpprintprec: 3;
  for i:0 thru 19 do (
    /* ldisplay( P_rodrigues[i] ), */
    /* ldisplay( P_recurrent[i] ), */
    /* ldisplay( P_generating[i] ), */
    /* ldisplay( P_hypergeometric[i] ), */
    /* printf(true, "~%----------------------~%") */
    results: addcol( results , [ measure_time(P_rodrigues[i]),
      measure_time(P_recurrent[i]),
      measure_time(P_generating[i]),
      measure_time(P_hypergeometric[i]) ] )
    );
texput(matrix, lambda([m], block(
    [rows : length(m), cols : length(m[1]), r, c, s: ""],
    s : sconcat("\\begin{pmatrix}", newline),
    for r:1 thru rows do (
       for c:1 thru cols do (
	  s : sconcat(s, tex1(m[r][c]), " & ")
       ),
       s : sconcat(s, "\\\\" , newline)
    ),
    sconcat(s, "\\end{pmatrix}" )
 )));
disp(sconcat( "\\[ ", tex1(results), " \\]"  ));

rodrigues__time=0.0006490000000000003 recurrent__time=0.0008319999999999994 generating_time=0.0006519999999999998 hypergeome_time=0.00048000000000000126 \[ \begin{pmatrix} {\it rodrigues} & 3.72 \times 10^{-4} & 2.88 \times 10^{-4} & 5.52 \times 10^{\ -4} & 8.49 \times 10^{-4} & 0.00105 & 0.0015 & 0.00191 & 0.00248 & 0.00326 & 0\ .00457 & 0.00478 & 0.00574 & 0.00683 & 0.00897 & 0.00913 & 0.0108 & 0.0149 & 0\ .0146 & 0.0155 & 0.0173 & \\ {\it recurrent} & 3.6 \times 10^{-5} & 3.7 \times 10^{-5} & 2.56 \times 10^{-4\ } & 5.36 \times 10^{-4} & 6.38 \times 10^{-4} & 7.88 \times 10^{-4} & 9.16 \ti\ mes 10^{-4} & 0.00108 & 0.00107 & 0.00146 & 0.00148 & 0.00166 & 0.00165 & 0.00\ 212 & 0.00187 & 0.00205 & 0.00219 & 0.0028 & 0.00242 & 0.00269 & \\ {\it generating} & 0.00607 & 0.00627 & 0.00589 & 0.00622 & 0.00603 & 0.0127 & \ 0.00626 & 0.00614 & 0.00637 & 0.00632 & 0.00618 & 0.00733 & 0.00634 & 0.00656 \ & 0.00654 & 0.00664 & 0.00706 & 0.00787 & 0.00695 & 0.00736 & \\ {\it hypergeometric} & 0.067 & 4.95 \times 10^{-4} & 8.56 \times 10^{-4} & 0.0\ 0131 & 0.00206 & 0.00306 & 0.00427 & 0.00653 & 0.0071 & 0.00874 & 0.0173 & 0.0\ 137 & 0.0161 & 0.0264 & 0.022 & 0.0316 & 0.0283 & 0.0383 & 0.0441 & 0.0422 & \\ \ \end{pmatrix} \]

Всё это выглядит более-менее пренебрежимо малым, однако, на первый взгляд, кажется, что гипергеометрический метод, а так же метод Родрига могут и подтормаживать.

2.2. Сделать трансляцию на язык Си программы на Maxima

Программа должна определять две функции, $f$ и $T$, являющиеся двумя разными представлениями одной и той же функции: символьным выражением и рядом Тейлора до степени 4

\[ f = \frac{1}{5+\sin(x)} \]

И сравнить результат от 0 до 1.0 с шагом 0.05.

load("gentran");
gentranlang: c;
f: 1/(5+ sin(x));
T(x) := part(taylor(f, x, 0,4), [1,2,3,4,5]);
print("#include<math.h>");
print("#include<stdio.h>");
print( "double f(double x){ return ");
gentran(eval(f));
print( "; }");
print( "");
print( "double T(double x){ return ");
gentran(eval(T(x)));
print("; }");
print("");
print("int main() { // return 0;");
print("for(double i=0; i<=1; i=i+0.05) { printf(\"%f\\t%f\\t%f\\t%f\\n\", i, f(i), T(i), f(i) - T(i)) ;}");
print("return 0;}");

#include<math.h> 
#include<stdio.h> 
double f(double x){ return  
1/(5+sin(x))
; } 

double T(double x){ return  
-22.0/9375.0*pow(x,4.0)+19.0/3750.0*pow(x,3.0)+pow(x,2.0)/125.0-1.0/25.0*x+1.0/
 5.0
; } 

int main() { // return 0; 
for(double i=0; i<=1; i=i+0.05) { printf("%f\t%f\t%f\t%f\n", i, f(i), T(i), f(\
i) - T(i)) ;} 
return 0;}

Okay, gentran should do it better than my crude solution here, but I am not sure how to do it better. Instead gentran_on does not seem to do what I expect it to do, and going through a stuff with maxima itself does not seem to make a lot of sense to me. It is nice to know that maxima supports , but until I have a piece of software to use this trick, I have no use for it.

2.3. Векторы состояния и линейные операторы

Написать программу, которая вычисляет результат действия произвольной функции операторов на произвольный вектор состояния.

Функция операторов уничтожения $a_{1},a_{2}$ , действующая на произвольный вектор состояния двумерного линейного гармонического осциллятора:

\[ |\psi> = \sum_{i=i}^{N}\sum_{j=i}^{M}C_{i,j}|i,j>, a_{1}|i,j>_{}= \sqrt{i}|i-1,j> , a_{2}|i,j> = \sqrt{j}|i,j-1> \]

Что-то у меня в формуле выше векторы бра и кет малость поломанные, ну да ладно.

Нужно посчитать $\sin(a_{1} + a_{2} * a_{1} )[w * |4,8> - (3/7)*|11,3>] $

Так, что тут можно сделать?

Что нужно вспомнить о квантовой механике человеку, который последний раз видел её больше 10 лет назад?

$|i>$ – это хитрое обозначение для вектора $(0\ 0\ 0\ \ldots\ 1\ \ldots\ ) $ в дискретном случае. Часто говорят, что это "система, в которой $i$ частиц", хотя на самом деле гораздо чаще это система, в которой единственная частица находится в $i$-раз возбуждённом состоянии.

Вектор состояния у меня будет sparse-матрицей в компьютерной смысле, вестимо.

Что такое $w$ в этой задаче? Окей, будем считать, что это символьный параметр, который должен в итоге войти в ответ.

В целом, хочется, например, получить результат в виде символьной матрицы?

input[4,8]  : w;
input[11,3] : -3/7;
display(arrayinfo(input));
/* Как задать input одним выражением? */
my_opfun[x]:= block([useless], 1);
a1(x) := block([retval] );
display(my_opfun[3]);

arrayinfo(input) = [hashed, 2, [4, 8], [11, 3]] my_opfun = 1 3

test1 : genmatrix( lambda( [n,m] , (if n>m then 0 else 1) * (n*m) ), 10, 10);
/*b : lreduce(cons, test1);*/
ldisp(test1);
ldisp(map(lambda([ROW],rest(ROW)), test1));
multiplies : genmatrix( lambda( [n,m], \sqrt(m)), 10,10);
ldisp(multiplies);
ldisp(matrixp(multiplies));
ldisp( test1 * multiplies );
ldisp( matrixp( test1 * multiplies) );
s1 : test1 * multiplies;
ldisp(matrixp(s1));
ldisp( addcol(submatrix(s1, 1), zeromatrix(length(s1),1)));

[ 1 2 3 4 5 6 7 8 9 10 ] [ ] [ 0 4 6 8 10 12 14 16 18 20 ] [ ] [ 0 0 9 12 15 18 21 24 27 30 ] [ ] [ 0 0 0 16 20 24 28 32 36 40 ] [ ] [ 0 0 0 0 25 30 35 40 45 50 ] (%t1) [ ] [ 0 0 0 0 0 36 42 48 54 60 ] [ ] [ 0 0 0 0 0 0 49 56 63 70 ] [ ] [ 0 0 0 0 0 0 0 64 72 80 ] [ ] [ 0 0 0 0 0 0 0 0 81 90 ] [ ] [ 0 0 0 0 0 0 0 0 0 100 ] [ 2 3 4 5 6 7 8 9 10 ] [ ] [ 4 6 8 10 12 14 16 18 20 ] [ ] [ 0 9 12 15 18 21 24 27 30 ] [ ] [ 0 0 16 20 24 28 32 36 40 ] [ ] [ 0 0 0 25 30 35 40 45 50 ] (%t2) [ ] [ 0 0 0 0 36 42 48 54 60 ] [ ] [ 0 0 0 0 0 49 56 63 70 ] [ ] [ 0 0 0 0 0 0 64 72 80 ] [ ] [ 0 0 0 0 0 0 0 81 90 ] [ ] [ 0 0 0 0 0 0 0 0 100 ] [ 3/2 ] [ 1 sqrt(2) sqrt(3) 2 sqrt(5) sqrt(6) sqrt(7) 2 3 sqrt(10) ] [ ] [ 3/2 ] [ 1 sqrt(2) sqrt(3) 2 sqrt(5) sqrt(6) sqrt(7) 2 3 sqrt(10) ] [ ] [ 3/2 ] [ 1 sqrt(2) sqrt(3) 2 sqrt(5) sqrt(6) sqrt(7) 2 3 sqrt(10) ] [ ] [ 3/2 ] [ 1 sqrt(2) sqrt(3) 2 sqrt(5) sqrt(6) sqrt(7) 2 3 sqrt(10) ] [ ] [ 3/2 ] [ 1 sqrt(2) sqrt(3) 2 sqrt(5) sqrt(6) sqrt(7) 2 3 sqrt(10) ] (%t3) [ ] [ 3/2 ] [ 1 sqrt(2) sqrt(3) 2 sqrt(5) sqrt(6) sqrt(7) 2 3 sqrt(10) ] [ ] [ 3/2 ] [ 1 sqrt(2) sqrt(3) 2 sqrt(5) sqrt(6) sqrt(7) 2 3 sqrt(10) ] [ ] [ 3/2 ] [ 1 sqrt(2) sqrt(3) 2 sqrt(5) sqrt(6) sqrt(7) 2 3 sqrt(10) ] [ ] [ 3/2 ] [ 1 sqrt(2) sqrt(3) 2 sqrt(5) sqrt(6) sqrt(7) 2 3 sqrt(10) ] [ ] [ 3/2 ] [ 1 sqrt(2) sqrt(3) 2 sqrt(5) sqrt(6) sqrt(7) 2 3 sqrt(10) ] (%t4) true [ 3/2 3/2 3/2 3/2 3/2 9/2 3/2 ] [ 1 2 3 8 5 6 7 2 27 10 ] [ ] [ 5/2 3/2 3/2 3/2 3/2 11/2 3/2 ] [ 0 2 2 3 16 2 5 2 6 2 7 2 54 2 10 ] [ ] [ 5/2 3/2 3/2 3/2 9/2 3/2 ] [ 0 0 3 24 3 5 3 6 3 7 3 2 81 3 10 ] [ ] [ 3/2 3/2 3/2 13/2 3/2 ] [ 0 0 0 32 4 5 4 6 4 7 2 108 4 10 ] [ ] [ 5/2 3/2 3/2 9/2 3/2 ] [ 0 0 0 0 5 5 6 5 7 5 2 135 5 10 ] (%t5) [ ] [ 5/2 3/2 11/2 3/2 ] [ 0 0 0 0 0 6 6 7 3 2 162 6 10 ] [ ] [ 5/2 9/2 3/2 ] [ 0 0 0 0 0 0 7 7 2 189 7 10 ] [ ] [ 15/2 3/2 ] [ 0 0 0 0 0 0 0 2 216 8 10 ] [ ] [ 3/2 ] [ 0 0 0 0 0 0 0 0 243 9 10 ] [ ] [ 5/2 ] [ 0 0 0 0 0 0 0 0 0 10 ] (%t6) true (%t7) true [ 3/2 3/2 3/2 3/2 3/2 9/2 3/2 ] [ 2 3 8 5 6 7 2 27 10 0 ] [ ] [ 5/2 3/2 3/2 3/2 3/2 11/2 3/2 ] [ 2 2 3 16 2 5 2 6 2 7 2 54 2 10 0 ] [ ] [ 5/2 3/2 3/2 3/2 9/2 3/2 ] [ 0 3 24 3 5 3 6 3 7 3 2 81 3 10 0 ] [ ] [ 3/2 3/2 3/2 13/2 3/2 ] [ 0 0 32 4 5 4 6 4 7 2 108 4 10 0 ] [ ] [ 5/2 3/2 3/2 9/2 3/2 ] [ 0 0 0 5 5 6 5 7 5 2 135 5 10 0 ] (%t8) [ ] [ 5/2 3/2 11/2 3/2 ] [ 0 0 0 0 6 6 7 3 2 162 6 10 0 ] [ ] [ 5/2 9/2 3/2 ] [ 0 0 0 0 0 7 7 2 189 7 10 0 ] [ ] [ 15/2 3/2 ] [ 0 0 0 0 0 0 2 216 8 10 0 ] [ ] [ 3/2 ] [ 0 0 0 0 0 0 0 243 9 10 0 ] [ ] [ 5/2 ] [ 0 0 0 0 0 0 0 0 10 0 ]

programmode: false;
  plot2d ([atan(x), erf(x), tanh(x)], [x, -5, 5], [y, -1.5, 2])$

a[2,3]: 1;
b[2,3]: 2;
s : ev(a+b);
ldisp(s);
programmode: false;
ldisp([1,2,3] + [5,5,5]);
map(ldisplay, a);

Loading home/lockywolf.maxima/maxima-init.mac (%t1) b + a (%t2) [6, 7, 8] map: improper argument: a – an error. To debug this try: debugmode(true);

programmode: false;
matchfix( "|" , ">" );
matchdeclare(a,numberp);
matchdeclare(b,numberp);
matchdeclare(c,true);
defmatch(ketp, c*|a,b>);
ldisp(ketp(|1,2>));
ldisp(ketp(alpha*|1,2>));
matchdeclare(a1,ketp);
matchdeclare(a2,ketp);
defmatch(qvectorp, a1 + a2);
ldisp(qvectorp(t*|1,2> + q*|1,2>));
tellsimp(a1 + a2, (op(a1)+op(a2)));
simp:true;
ldisp( ev(t*|1,2> + q*|1,2>));

Loading home/lockywolf.maxima/maxima-init.mac read and interpret /tmp/babel-cTZEIJ/maxima-VXgAxR.max set_tex_environment_default("\\[ "," \\]") [\[ , \]] programmode:false false matchfix("|",">")

matchdeclare(a,numberp) done matchdeclare(b,numberp) done matchdeclare(c,true) done defmatch(ketp,c*|a,b>) ketp ldisp(ketp(|1,2>)) (%t9) [c = 1, b = 2, a = 1] [%t9] ldisp(ketp(alpha*|1,2>)) (%t10) [c = alpha, b = 2, a = 1] [%t10] matchdeclare(a1,ketp) done matchdeclare(a2,ketp) done defmatch(qvectorp,a1+a2) qvectorp ldisp(qvectorp(t*|1,2>+q*|1,2>)) (%t14) false [%t14] tellsimp(a1+a2,op(a1)+op(a2)) tellsimp: warning: rule will treat '+' as noncommutative and nonassociative. [+rule1, simplus] simp:true true ldisp(ev(t*|1,2>+q*|1,2>)) (%t17) |1, 2> t + |1, 2> q [%t17] gnuplot_close() quit()

load(partition);
matchdeclare (pp, partition_expression("+",constantp,0,"+",g,'ANSpp));
tellsimp(foo(pp),ANSpp);
declare([a,b,c],constant);
disp([foo(a+b+c+x+y+3),foo(w),foo(b), foo(a*b+x),foo(a*b)]);

Loading home/lockywolf.maxima/maxima-init.mac [g(c + b + a + 3, y + x), foo(w), foo(b), g(a b, x), foo(a b)]

programmode: false;
nolabels: false;
matchdeclare([ann,bnn], lambda([r],atom(r) and not numberp(r)),
	    cna,lambda([r], not atom(r)),
	    numonly,numberp);
defmatch(m1,h(ann,bnn,numonly) );
defrule(r1,h(ann,bnn,numonly),['numonlyA = numonly,'bnnB = bnn,'annC = ann])$
ldisp(m1(h(q,r,34)));
ldisp(r1(h(q,r,34)));
dadfa;
:lisp (symbol-plist '$ann)

Loading home/lockywolf.maxima/maxima-init.mac read and interpret /tmp/babel-cTZEIJ/maxima-aoIKhA.max set_tex_environment_default("\\[ "," \\]") [\[ , \]] programmode:false false nolabels:false false matchdeclare([ann,bnn],lambda([r],atom(r) and not numberp(r)),cna, lambda([r],not atom(r)),numonly,numberp) done defmatch(m1,h(ann,bnn,numonly)) m1 defrule(r1,h(ann,bnn,numonly),['numonlyA = numonly,'bnnB = bnn,'annC = ann]) ldisp(m1(h(q,r,34))) (%t8) [numonly = 34, bnn = r, ann = q] [%t8] ldisp(r1(h(q,r,34))) (%t9) [numonlyA = 34, bnnB = r, annC = q] [%t9] dadfa dadfa (MPROPS (NIL MATCHDECLARE (((LAMBDA (0 /tmp/babel-cTZEIJ/maxima-aoIKhA.max SRC)) ((MLIST (0 /tmp/babel-cTZEIJ/maxima-aoIKhA.max SRC)) $R) ((MAND (0 /tmp/babel-cTZEIJ/maxima-aoIKhA.max SRC)) (($ATOM (0 /tmp/babel-cTZEIJ/maxima-aoIKhA.max SRC)) $R) ((MNOT (0 /tmp/babel-cTZEIJ/maxima-aoIKhA.max SRC)) (($NUMBERP (0 /tmp/babel-cTZEIJ/maxima-aoIKhA.max SRC)) $R))))))) gnuplot_close() quit()

After a bit of searching, I found this code, of Leo Butler 2010. (Michael Talon telling me at Maxima's mailing list.) It should be not hard to tweak it for my purposes?

  tellsimp(R.L,1+L.R);
  tellsimp(R^^i . L, 1+(R^^(i-1) . L) . R);
  declare(h,integer);
  declare(h,scalar);
  matchdeclare(m, lambda([t1],featurep(t1,integer)),
	       n, lambda([t2],featurep(t2,integer)),
	       u, scalarp);
  tellsimp(L.ket(n), sqrt(n)*ket(n-1));
  tellsimp(L.(u*ket(n)), u*sqrt(n)*ket(n-1));
  tellsimp(L^^m . (u*ket(n)), u*sqrt(n)*(L^^(m-1) . ket(n-1)));
  tellsimp(R.ket(n), sqrt(n+1)*ket(n+1));
  tellsimp(R.(u*ket(n)), u*sqrt(n+1)*ket(n+1));
  tellsimp(R^^m . (u*ket(n)), u*sqrt(n+1)*(R^^(m-1) . ket(n+1)));
  L.ket(h);
  L.ket(2);
  L.L.ket(h);
  L^^3 . ket(5);
ket(5) + ket(5);
R.ket(5) ;
ket(1,2) + ket(1,2);
L.ket(0);
(L + R) . ket(5);
telsimp(R+L, L + R);
L + R;
load("functs");
matchdeclare( qoperator, nonzeroandfreeof(ket) );
defmatch(qoperatorp, qoperator);
matchdeclare( qvector, lambda([t3], not(freeof(ket, t3))) );
matchdeclare( p, ?mplusp );
tellsimp( (p) . qvector,
	   first(p).qvector + rest(p).qvector);
(L + R) . ket(5);
(8*L + 7*R) . ket(5);
sin(8*L + 7*R) . ket(5);
matchdeclare( ker, lambda([t4], atom(t4) and not(numberp(t4))));
/*telsimp(testfun, ker(qoperator));*/
/*  defrule(r, ker(qoperator), powerseries(ker(qoperator), qoperator, 0));*/
matchdeclare(oppower, lambda([t5], integerp(t5) and is(t5>0)));
tellsimp( L^oppower, L^^oppower);
tellsimp( R^oppower, R^^oppower);
(R^5 . L^5).ket(5);
sin(8*L + L) . ket(5);
simpsum: true;
tellsimp( ker(qoperator) . qvector,
	  powerseries(ker(qoperator), qoperator, 0) . qvector) ;
simpsum: false;
mysum : op(sum(i^2,i,0,inf));
ldisplay(mysum);
matchdeclare(sump, lambda([t6], not(atom(t6)) and is(equal(op(t6),mysum))));
defmatch(t7,sump);
t7(sum(i,i,0,inf));
defmatch(t8, sump . qvector );
tellsimp( sump . qvector, apply(sum, cons(ratexpand(first(sump)) . qvector, rest(args(sump)))));
simpsum: false;
sin(8*L + L) . ket(5);
t8( sum(L^i,i,0,inf) . ket(5) );
(L).ket(5);

Loading home/lockywolf.maxima/maxima-init.mac read and interpret /tmp/babel-sPWpfR/maxima-NOwpwS.max set_tex_environment_default("\\[ "," \\]") [\[ , \]] tellsimp(R . L,1+L . R) [.rule1, simpnct] tellsimp(R^^i . L,1+(R^^(i-1) . L) . R) [.rule2, .rule1, simpnct] declare(h,integer) done declare(h,scalar) done matchdeclare(m,lambda([t1],featurep(t1,integer)),n, lambda([t2],featurep(t2,integer)),u,scalarp) done tellsimp(L . ket(n),sqrt(n)*ket(n-1)) [.rule3, .rule2, .rule1, simpnct] tellsimp(L . (u*ket(n)),u*sqrt(n)*ket(n-1)) [.rule4, .rule3, .rule2, .rule1, simpnct] tellsimp(L^^m . (u*ket(n)),u*sqrt(n)*L^^(m-1) . ket(n-1)) [.rule5, .rule4, .rule3, .rule2, .rule1, simpnct] tellsimp(R . ket(n),sqrt(n+1)*ket(n+1)) [.rule6, .rule5, .rule4, .rule3, .rule2, .rule1, simpnct] tellsimp(R . (u*ket(n)),u*sqrt(n+1)*ket(n+1)) [.rule7, .rule6, .rule5, .rule4, .rule3, .rule2, .rule1, simpnct] tellsimp(R^^m . (u*ket(n)),u*sqrt(n+1)*R^^(m-1) . ket(n+1)) [.rule8, .rule7, .rule6, .rule5, .rule4, .rule3, .rule2, .rule1, simpnct] L . ket(h) ket(h - 1) sqrt(h) L . ket(2) sqrt(2) ket(1) L . L . ket(h) ket(h - 2) sqrt(h - 1) sqrt(h) L^^3 . ket(5) 2 sqrt(3) sqrt(5) ket(2) ket(5)+ket(5) 2 ket(5) R . ket(5) sqrt(6) ket(6) ket(1,2)+ket(1,2) 2 ket(1, 2) L . ket(0) 0 (L+R) . ket(5) (R + L) . ket(5) telsimp(R+L,L+R) telsimp(R + L, R + L) L+R R + L load("functs") /usr/share/maxima/5.47.0/share/simplification/functs.mac matchdeclare(qoperator,nonzeroandfreeof(ket)) done defmatch(qoperatorp,qoperator) defmatch: evaluation of atomic pattern yields: qoperator qoperatorp matchdeclare(qvector,lambda([t3],not freeof(ket,t3))) done matchdeclare(p,mplusp) done tellsimp(p . qvector,first(p) . qvector+rest(p) . qvector) [.rule9, .rule8, .rule7, .rule6, .rule5, .rule4, .rule3, .rule2, .rule1, simpnct] (L+R) . ket(5) sqrt(6) ket(6) + sqrt(5) ket(4) (8*L+7*R) . ket(5) 7 sqrt(6) ket(6) + 8 sqrt(5) ket(4) sin(8*L+7*R) . ket(5) sin(7 R + 8 L) . ket(5) matchdeclare(ker,lambda([t4],atom(t4) and not numberp(t4))) done matchdeclare(oppower,lambda([t5],integerp(t5) and is(t5 > 0))) done tellsimp(L^oppower,L^^oppower) [^rule1, simpexpt] tellsimp(R^oppower,R^^oppower) [^rule2, ^rule1, simpexpt] (R^5 . L^5) . ket(5) 120 ket(5) sin(8*L+L) . ket(5) sin(9 L) . ket(5) simpsum:true true tellsimp(ker(qoperator) . qvector, powerseries(ker(qoperator),qoperator,0) . qvector) [.rule10, .rule9, .rule8, .rule7, .rule6, .rule5, .rule4, .rule3, .rule2, .rule1, simpnct] simpsum:false false mysum:op(sum(i^2,i,0,inf)) sum ldisplay(mysum) (%t44) mysum = sum [%t44] matchdeclare(sump,lambda([t6],not atom(t6) and is(equal(op(t6),mysum)))) done defmatch(t7,sump) defmatch: evaluation of atomic pattern yields: sump t7 t7(sum(i,i,0,inf)) inf == \ [sump = > i] / == i = 0 defmatch(t8,sump . qvector) t8 tellsimp(sump . qvector, apply(sum,cons(ratexpand(first(sump)) . qvector,rest(args(sump))))) [.rule11, .rule10, .rule9, .rule8, .rule7, .rule6, .rule5, .rule4, .rule3, .rule2, .rule1, simpnct] simpsum:false false sin(8*L+L) . ket(5) inf == 2 i1 + 1 i1 2 i1 + 1 \ L (- 1) 9 > (----------------------–—) . ket(5) / (2 i1 + 1)! == i1 = 0 t8(sum(L^i,i,0,inf) . ket(5)) false L . ket(5) sqrt(5) ket(4) gnuplot_close() quit()

matchdeclare(m, lambda([t],featurep(t,integer)), /* power */
  n1, lambda([t],featurep(t,integer)), /* vector1 */
  n2, lambda([t],featurep(t,integer)), /* vector2 */
  u, scalarp); /* argument */
tellsimp(a1.ket(n1,n2), sqrt(n1)*ket(n1-1,n2));
tellsimp(a1.(u*ket(n1,n2)), u*sqrt(n1)*ket(n1-1,n2));
tellsimp(a1^^m . (u*ket(n1,n2)), u*sqrt(n1)*(a1^^(m-1) . ket(n1-1,n2)));
/* tellsimp(a1^^m . (u*ket(n1,n2)), a1^^(m-1). a1 . ket(n1-1,n2)); */

a1.ket(2,3);
a1^^3 . ket(5,1);

tellsimp(a2.ket(n1,n2), sqrt(n2)*ket(n1,n2-1));
tellsimp(a2.(u*ket(n1,n2)), u*sqrt(n2)*ket(n1,n2-1));
tellsimp(a2^^m . (u*ket(n1,n2)), u*sqrt(n2)*(a2^^(m-1) . ket(n1,n2-1)));

a2.ket(2,3);
a2^^3 . ket(1,5);

a2.a1.(w * ket(4,8) - (3/7) * ket(11,3));
load("diag");
mat_function(sin,(a1 + a2*a1)(w * ket(4,8) - (3/7) * ket(11,3)));
A : matrix([2,4],[1,2]);
M : mat_function(sin,t*A);

Loading home/lockywolf.maxima/maxima-init.mac read and interpret /tmp/babel-cTZEIJ/maxima-3v9KzK.max set_tex_environment_default("\\[ "," \\]") [\[ , \]] matchdeclare(m,lambda([t],featurep(t,integer)),n1, lambda([t],featurep(t,integer)),n2, lambda([t],featurep(t,integer)),u,scalarp) done tellsimp(a1 . ket(n1,n2),sqrt(n1)*ket(n1-1,n2)) [.rule1, simpnct] tellsimp(a1 . (u*ket(n1,n2)),u*sqrt(n1)*ket(n1-1,n2)) [.rule2, .rule1, simpnct] tellsimp(a1^^m . (u*ket(n1,n2)),a1^^(m-1) . a1 . ket(n1-1,n2)) [.rule3, .rule2, .rule1, simpnct] a1 . ket(2,3) sqrt(2) ket(1, 3) a1^^3 . ket(5,1) Maxima encountered a Lisp error: Binding stack exhausted. PROCEED WITH CAUTION. Automatically continuing. To enable the Lisp debugger set debugger-hook to nil. tellsimp(a2 . ket(n1,n2),sqrt(n2)*ket(n1,n2-1)) [.rule4, .rule3, .rule2, .rule1, simpnct] tellsimp(a2 . (u*ket(n1,n2)),u*sqrt(n2)*ket(n1,n2-1)) [.rule5, .rule4, .rule3, .rule2, .rule1, simpnct] tellsimp(a2^^m . (u*ket(n1,n2)),u*sqrt(n2)*a2^^(m-1) . ket(n1,n2-1)) [.rule6, .rule5, .rule4, .rule3, .rule2, .rule1, simpnct] a2 . ket(2,3) sqrt(3) ket(2, 2) a2^^3 . ket(1,5) 2 sqrt(3) sqrt(5) ket(1, 2) a2 . a1 . (w*ket(4,8)-(3/7)*ket(11,3)) 3 ket(11, 3) a2 . a1 . (ket(4, 8) w - -------–—) 7 load("diag") /usr/share/maxima/5.47.0/share/contrib/diag.mac mat_function(sin,(a1+a2*a1)(w*ket(4,8)-(3/7)*ket(11,3))) 3 ket(11, 3) mat_function(sin, a1 a2 + a1(ket(4, 8) w - -------–—)) 7 A:matrix([2,4],[1,2]) [ 2 4 ] [ ] [ 1 2 ] M:mat_function(sin,t*A) [ sin(4 t) ] [ ---–— sin(4 t) ] [ 2 ] [ ] [ sin(4 t) sin(4 t) ] [ ---–— ---–— ] [ 4 2 ] gnuplot_close() quit()

Что-то я возился-возился, ни черта не понял. Задачка должна быть простая относительно. Но что-то я ухожу всё дальше и дальше в дебри.

Если реализовать операторы как "переменные", то получается ерунда с разложением их в ряд по степеням, ведь там нужна матричная степень, а не "обычная". Если их реализовывать как кастомные операторы, то нельзя сотворить выражение, не содержащее аргументов. mat_function отказывается раскладывать выражение в ряд, если в нём нет матриц, но с матрицами плохо, потому что кет-вектор плохо представляется в качестве одномерного вектора.

И это простая задача! Ничего не понимаю. Если в последнем случае поменять simpsum:false на simpsum:true, система впадает в бесконечный цикл, что плохо, :(. Не знаю, что на этом дальше делать.

2.4. Векторы состояния и линейные операторы

Написать "алгебру операторов гармонического осциллятора", для приведения выражений с операторами к каноническому виду.

Тестовый пример:

-(76/3m)LRRLLR + 23LR + … +5RR

Все операторы R должны стоять слева от операторов L.

Ну здесь-то можно прямолинейно применить "алгебру гармонического осциллятора" из предыдущего упражнения?

tellsimpafter(R.L,I+L.R);
tellsimpafter(R^^i . L, I+(R^^(i-1) . L) . R);
tellsimpafter(R.I, R);
tellsimpafter(I.R, R);
tellsimpafter(L.I, L);
tellsimpafter(I.L, L);

a: ev((-(76/(3*m))*L.R.R.L.L.R + 23*L.R + L.L.L.R.L.R.L + 5*R.R), expand);
expand(a);

Loading home/lockywolf.maxima/maxima-init.mac read and interpret /tmp/babel-sPWpfR/maxima-yHv61R.max set_tex_environment_default("\\[ "," \\]") [\[ , \]] tellsimpafter(R . L,I+L . R) [.rule1, simpnct] tellsimpafter(R^^i . L,I+(R^^(i-1) . L) . R) [.rule2, .rule1, simpnct] tellsimpafter(R . I,R) [.rule3, .rule2, .rule1, simpnct] tellsimpafter(I . R,R) [.rule4, .rule3, .rule2, .rule1, simpnct] tellsimpafter(L . I,L) [.rule5, .rule4, .rule3, .rule2, .rule1, simpnct] tellsimpafter(I . L,L) [.rule6, .rule5, .rule4, .rule3, .rule2, .rule1, simpnct] a:ev(-(76/(3*m))*L . R . R . L . L . R+23*L . R+L . L . L . R . L . R . L +5*R . R,expand) <2> <2> 76 (L . R . L . R) <2> <4>

-------------------–— + 5 R + L . R + 23 (L . R) 3 m <2> <2> <3>
- L . (L . R) + L

expand(a) <2> <2> 76 (L . R . L . R) <2> <4>

-------------------–— + 5 R + L . R + 23 (L . R) 3 m <2> <2> <3>
- L . (L . R) + L

gnuplot_close() quit()

Опять чудеса. Ладно, спросим в мейлинг листе Максимы, что вообще происходит.

2.5. Большая задача. Получить спектральное уравнение для малых аксиальных возмущений решения Шварцшильда в ОТО.

В общем, пока я не прочитаю хотя бы введение в ОТО, маловероятно, что я решу эту задачу. На выбор есть ещё четыре задачи, у которых я не знаю всех терминов, входящих в условия.

3. Заключение

В общем, проект зафейлился. Может, лет через 5 я выучу достаточно, чтобы понимать, о чём речь, и вернусь к нему.

1. TODO todos [0/1]