Fri, Apr. 24th, 2009, 01:15 pm
Что с вашей Максимой?

Обнаружил сегодня вот такой забавный баг в Максиме:
Пусть имеются две функции, например,
f1:4*x^3-4*y;
f2:4*y^3-4*x;
Далее, пытаюсь её решить
solve([f1=0,f2=0],[x,y]);
Как и полагается, она выводит 9 вариантов решений, из которых 6 - комплексные. Далее, я хочу комплексные решения убрать:
solve([f1=0,f2=0],[x,y]),realonly:true;
она убирает все те комплексные варианты, которые в явном виде содержат i. Однако те, которые содержат (-1)^(1/4) почему-то остаются.

Это неправильно же, мы все знаем, что (-1)^(1/4)=sqrt(i)=sqrt(2)/2+(sqrt(2)/2)i, то есть ниразу не действительное число. Как бы это Максиме объяснить?
 Leave a comment | 17 comments

Fri, Apr. 24th, 2009, 05:03 am
[info]beshenov

Это не баг. В документации явно сказано, что realonly исключает решения с %i. Там не сказано ничего про комплексные решения.

Я бы сделал так (на скорую руку):

filter(pred, lst) :=
  if lst=[] then [] else
  block ([x : first(lst), xs : rest(lst)],
    if pred(x) then cons(x,filter(pred,xs)) else filter(pred,xs)
  )$

is_real(x) := is(trigsimp(imagpart(x))=0)$

solve_real([args]) :=
    block (
      [sol : apply('solve, args)],

      if sol#[] then
        if listp(first(sol)) then
          delete ([],
                  map (lambda([xs],
                              filter(lambda([x],is_real(rhs(x))), xs)), 
                       sol))
        else
          filter(lambda([x],is_real(rhs(x))), sol)
    )$


Результат:

solve_real([f1=0,f2=0],[x,y]) =>
    [[x = - 1, y = - 1], [x = 1, y = 1], [x = 0, y = 0]]


trigsimp(imagpart(x))=0 должно обезопасить от casus irreducibilis, например.
Thread - Reply

Fri, Apr. 24th, 2009, 05:34 am
[info]ash_rabbi

Ну блин... Это ж сложно так, и в imath не вставить... То есть, вставить можно, и даже работает, но output вставляется не тот. Может быть, есть какой-нибудь ключ, чтобы она все комплексные числа пыталась приводить к виду a+bi хотя бы?
Parent - Thread - Reply

Fri, Apr. 24th, 2009, 08:40 am
[info]beshenov

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

Есть функция rectform. Я привел просто общее решение проблемы, которое можно записать в foo.mac и использовать через load("foo.mac").

На практике я бы просто написал после solve([f1=0,f2=0],[x,y]); строчку map ('rectform, %);
Parent - Thread - Reply

Fri, Apr. 24th, 2009, 08:43 am
[info]beshenov

Или сразу solve([f1=0,f2=0],[x,y]), rectform;
Parent - Thread - Reply

Fri, Apr. 24th, 2009, 10:36 am
[info]ash_rabbi

Непосредственно в таком виде оно даёт просто все корни в нормальной форме. Так же, как и rectform(solve([f_x=0,f_y=0],[x,y])); Я попытался сделать rectform(solve([f_x=0,f_y=0],[x,y])),realonly:true; но почему-то в таком виде комплексные корни не убираются, несмотря на то, что там вполне себе есть i.
Parent - Thread - Reply

Fri, Apr. 24th, 2009, 03:03 pm
[info]beshenov

Так флаг realonly действует на работу solve (algsys, точнее), а rectform --- при упрощении выражений в том виде, в каком их уже выдала функция algsys.
Parent - Thread - Reply

Sat, Apr. 25th, 2009, 01:49 am
[info]ash_rabbi

Блин, а нет какой-нибудь функции, которая могла бы просто отрезать из вывода все корни с i? Просто хотелось бы это сделать в одну строчку, потому что иначе imath не работает по-человечески...
Parent - Reply

Sat, Apr. 25th, 2009, 04:15 am
[info]ash_rabbi

Короче, на этот раз выкрутился так: realp(x):=is(imagpart(x)=[0=0,0=0])$sublist(solve([f_x=0,f_y=0],[x,y]),realp); Это как бы две строчки, но imath не заметил подмены.

Жуткий костыль, конечно, но в данной конкретной задаче справился. А на будущее надо, я думаю, придумать нормальный предикат и внедрить его в максиму. Потому что, к примеру, oddp, evenp, и ratnump там есть, почему нет realp?
Parent - Thread - Reply

Sat, Apr. 25th, 2009, 04:53 am
[info]beshenov

Я уже практически всё написал. Особенности и ограничения отдельного неофициального интерфейса мне не интересны. Тем более что есть load() и пользовательские файлы инициализации.

Просто imagpart(x)=0 не годится (и тем более freeof (%i, ...)). Возьмем какой-нибудь неприводимый полином с вещественными корнями, допустим 3 x^3 - 3 x + 1.

(%i1) solve (3*x^3  - 3*x + 1);
(%o1) [x = (sqrt(3)*%i/2-1/2)/(3*(3^-(3/2)*%i/2-1/6)^(1/3))
         +(3^-(3/2)*%i/2-1/6)^(1/3)*(-sqrt(3)*%i/2-1/2),
       x = (3^-(3/2)*%i/2-1/6)^(1/3)*(sqrt(3)*%i/2-1/2)
         +(-sqrt(3)*%i/2-1/2)/(3*(3^-(3/2)*%i/2-1/6)^(1/3)),
       x = (3^-(3/2)*%i/2-1/6)^(1/3)+1/(3*(3^-(3/2)*%i/2-1/6)^(1/3))]

(%i2) map (lambda ([x], is(imagpart(rhs(x))=0)), %o1);
(%o2) [false,false,true]

(%i3) map (lambda ([x], is(trigsimp(imagpart(rhs(x)))=0)), %o1);
(%o3) [true,true,true]
Parent - Thread - Reply

Sun, Apr. 26th, 2009, 08:54 am
[info]ash_rabbi

Ну так может быть имеет смысл сделать в виде патча и послать разработчикам, чтобы они внедрили непосредственно в дефолтную сборку? Тогда оно будет работать везде, во всех интерфейсах. Кроме того, realonly в текущем виде, по-моему, где-то настолько же полезен, как тот томагавк, который может вскрывать черепа только диким цесаркам не младше одного, но не старше двух лет. По крайней мере, мне трудно представить такую задачу, в которой было бы необходимо удалить все решения с %i, но при этом оставить решения с (-1)^(1/4) и (-1)^(3/4).
Parent - Thread - Reply

Sun, Apr. 26th, 2009, 09:41 am
[info]beshenov

Я спрошу на листе и добавлю, если подобные изменения будут одобрены.
Parent - Reply

Fri, Apr. 24th, 2009, 08:58 am
[info]beshenov

Уточнение:

is_real(x) := is(trigsimp(imagpart(x))=0)$

solve_real([args]) :=
    block (
      [sol : apply('solve, args)],

      if sol#[] then
        if listp(first(sol)) then
          delete ([],
                  map (lambda([xs],
                              sublist(xs, lambda([x],is_real(rhs(x))))), 
                       sol))
        else
          filter(lambda([x],is_real(rhs(x))), sol)
      else []
    )$

f1 : 4*x^3 - 4*y$
f2 : 4*y^3 - 4*x$

solve([f1=0,f2=0], [x,y]) =>
    [[x = - 1, y = - 1], [x = 1, y = 1], [x = 0, y = 0]]
Parent - Reply

Fri, Apr. 24th, 2009, 08:59 am
[info]beshenov

(Перепутал, "filter" уже есть и называется sublist.)
Parent - Reply

Sun, May. 3rd, 2009, 03:35 am
[info]laplandian.livejournal.com

Я бы так написал:

solve([f1=0,f2=0],[x,y]),rectform;sublist(%,lambda([e],freeof(%i,e)));
Thread - Reply

Sun, May. 3rd, 2009, 07:36 am
[info]ash_rabbi

В твоём варианте тоже какая-то чушь получается с неприводимыми уравнениями. Например, пример тов. [info]beshenov:
ff:3*x^3- 3*x + 1;
solve(ff);
trigsimp(%);

даёт следующее:

[x=-((%i-sqrt(3))^(2/3)*(3*sqrt(3)*%i+3)-3^(5/6)*6^(2/3)*%i+3^(1/3)*6^(2/3))/(6*3^(1/6)*6^(1/3)*(%i-sqrt(3))^(1/3)),x=((%i-sqrt(3))^(2/3)*(3*sqrt(3)*%i-3)-3^(5/6)*6^(2/3)*%i-3^(1/3)*6^(2/3))/(6*3^(1/6)*6^(1/3)*(%i-sqrt(3))^(1/3)),x=(3*(%i-sqrt(3))^(2/3)+3^(1/3)*6^(2/3))/(3*3^(1/6)*6^(1/3)*(%i-sqrt(3))^(1/3))]

То есть, мы видим, что действительных корней у него нет. Тем не менее, твой вариант
solve([ff=0],[x]),rectform$sublist(%,lambda([e],freeof(%i,e)));
даёт
[x=2*3^(-1/2)*cos((%pi-atan(3*3^(-3/2)))/3)]
притом он не упрощается:
trigsimp(%)
[x=(2*cos((5*%pi)/18))/sqrt(3)]
если его подставить в исходное уравнение
ff,x=(2*cos((5*%pi)/18))/sqrt(3);
получится
(8*cos((5*%pi)/18)^3)/sqrt(3)-(6*cos((5*%pi)/18))/sqrt(3)+1
то есть, это совершенно даже не корень уравнения и вообще непонятно, откуда он взялся.

Надо как-то привинтить is(imagpart(trigsimp(e))=0), но я не достаточно хорошо знаю лисп, чтобы понять, как. Если применить imagpart прямо к сублисту, получится не 0, а 0=0,0=0,0=0 итп столько раз, сколько было переменных. Как-нибудь можно применить фильтр только к тому, что находится между знаком равенства и запятой?
Parent - Thread - Reply

Sun, May. 3rd, 2009, 10:17 am
[info]beshenov

Есть lhs / rhs, которые выделяют части выражений вида a = b.

Lisp НЕ НАДО знать. Надо написать всё на встроенном языке Maxima и записать в свой файл инициализации в домашней директории.

Я писал о проблеме в список рассылки. Надо бы багрепорт, чтобы не потерялось. (По пути выяснилось, что для корней неприводимых полиномов abs выдает разные результаты в зависимости от того, приведены ли они явно к rectform. А это уже явно баг.)
Parent - Thread - Reply

Sun, May. 3rd, 2009, 12:24 pm
[info]ash_rabbi

Ну, насколько я правильно понял, максима может чуть менее чем всё, что может лисп. А вот в свой личный файл инициализации я ничего прятать не могу, потому что здешние преподы ничего кроме матлаба в глаза не видели, а мне совершенно не интересно им объяснять, с чего это у меня такая-то команда работает, а у них, когда они пытаются воспроизвести мой код - нет.

Про rhs я думал, конечно. Придумал вот такое:
solve([f]),rectform$ sublist(trigsimp(%),lambda([e],is(imagpart(rhs(e))=0)));
Но это не работает, когда более одной переменной.

Лучше бы, конечно, realonly в самой максиме поправили. Но я не знаю, что там за сообщество и как быстро оно реагирует...
Parent - Reply