Метод Кардано и Mathcad

[DEP]

Использовал Mathcad 15

Кубическое уравнение
$a_3x^3 + a_2x^2 + a_1x +a_0 = 0$ (1)
Делаю подстановку $x = y - \frac{a_2}{3a_3}$ получаю уравнение
$y^3 + qy + p = 0$ (2)
- это уравнение подстановкой $y = t - \frac{p}{3t}$ привожу к виду
$t^6 + qt^3 + \frac{p^3}{27}$ (3)
где $t^3 = w$
И затем решаю уравнение $w^2 + qw + \frac{p^3}{27}$ (4)
в Mathcad
1. по формуле $w_{1,2} = \frac{- q}{2} \pm \sqrt{\frac{q^2}{4} + \frac{p^3}{27}}$
2. с помощью оператора solve

Оба результата совпадают, корни $w_{1,2}$ комплексно сопряжённые, подставляя их в (4) получаю ноль. Далее в Mathcad по формуле Муавра извлекаю кубические корни из $w_1$ и $w_2$, получаю два вектора $T1$ и $T2$, вектор $T1$ содержит по три корня уравнения (3) комплексно сопряжённые к корням вектора $T2$.
Так как подстановка $y = t - \frac{p}{3t}$ приводится к виду $y = \sqrt[3]{\frac{- q}{2} + \sqrt{\frac{q^2}{4} + \frac{p^3}{27}}} + \sqrt[3]{\frac{- q}{2} - \sqrt{\frac{q^2}{4} + \frac{p^3}{27}}}$, это значит что элементы из векторов $T1$ и $T2$ я должен, комбинируя друг с другом и складывая, подставить в уравнение (2) и узнать какой из корней обнулит его. Так и делаю, из 9 сумм $y_n = T1_i + T2_j$ только одна является корнем уравнения (2), вышло что суммирую два комплексно сопряжённых корня и получаю один действительный.

Этот один действительный корень $y_1$ уравнения (2) я с помощью подстановки $x_1 = y_1 - \frac{a_2}{3a_3}$ подставляю в уравнение (1) и.... оно не равно нулю... а равно ненулевому остатку от деления уравнения (1) по схеме Горнера на двучлен $x - x_1$. Почему ?

И все 9 выше указанных псевдо-корней $x_n = y_n - \frac{a_2}{3a_3}$ тоже не обнуляют уравнение (1)... Точность вычислений поставил на максимум (307 знаков), значит ни один из коэффициентов не обнуляется... Что же не так делаю, где ошибка ?

[Pphantom]

DEP в сообщении #1240270 писал(а):

Точность вычислений поставил на максимум (307 знаков)

Так Вы решаете уравнение в общем виде или с какими-то конкретными значениями коэффициентов? Если второе - приведите их, найденные корни и значения многочлена в этих точках.

[DEP]

Pphantom в сообщении #1240279 писал(а):

Так Вы решаете уравнение в общем виде или с какими-то конкретными значениями коэффициентов? Если второе - приведите их, найденные корни и значения многочлена в этих точках.

Пытаюсь посчитать колебательный переходный процесс в RLC цепи 3-го порядка.

$a_3 = 3.48\cdot10^{-18} , a_2 = 1.76\cdot10^{-12} , a_1 = -1.54\cdot10^{-6} , a_0 = 0.09$
$p = -5.278\cdot10^{11} , q = 1.100468425\cdot10^{17}$

Кажется разобрался, наверное из-за очень больших и коэфициентов p и q внутри Mathcad происходит странное округление промежуточных результатов, и подстановка корня в уравнение даёт не ноль а что-то вида 736+170.660833i и это меня сбивало с толку.

Получилось 3 действительных корня, как и через оператор solve, но отличаются в каком-то десятом знаке, и вот...

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

[svv]

DEP в сообщении #1240299 писал(а):

$a_3 = 3.48\cdot10^{-18} , a_2 = 1.76\cdot10^{-12} , a_1 = -1.54\cdot10^{-6} , a_0 = 0.09$
$p = -5.278\cdot10^{11} , q = 1.100468425\cdot10^{17}$

Кажется разобрался, наверное из-за очень больших коэфициентов $p$ и $q$ внутри Mathcad происходит странное округление промежуточных результатов

В данном случае напрашивается подстановка $x=10^6 z$. Уравнение для $z$ будет иметь вид
$3.48z^3+1.76z^2-1.54z+0.09=0$
Совсем другое дело, правда?

[Pphantom]

DEP в сообщении #1240299 писал(а):

Кажется разобрался, наверное из-за очень больших и коэфициентов p и q внутри Mathcad происходит странное округление промежуточных результатов

Именно. Рецепт борьбы уже предложил svv.

[GAA]

DEP в сообщении #1240299 писал(а):

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

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

Mathcad 13 (символьный движок под Maple) ищет решения данного уравнения корректно (см. картинку)
Mathcad 15 (символьный движок MuPAD) сейчас нет под рукой.

-- Сб 12.08.2017 20:00:39 --

Найденные в Maple 7.0 корни (Digits=32)

-982797.21044809201225727922804280, 63655.263354112093488703161062167, 413394.82065719830957317376813005

[DEP]

svv в сообщении #1240305 писал(а):

В данном случае напрашивается подстановка $x=10^6 z$. Уравнение для $z$ будет иметь вид
$3.48z^3+1.76z^2-1.54z+0.09=0$
Совсем другое дело, правда?

Спасибо, проверю.

GAA в сообщении #1240308 писал(а):

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

Вобщем разобрался, чтобы у кубического уравнения было 2 комплексно сопряжённых корня - дискриминант уравнения (4) должен быть >0, тоесть по Формуле Муавра кубические корни из действительных $w_{1,2}$ не будут взаимно комплексно сопряжёнными, и комлексность корней $y_n = T1_i + T2_j$ не уничтожится при сложении.
Решая уравнение в котором один параметр схемы переменный, я нахожу его значение при котором дискриминант равен нулю, и варьируя этим параметром легко получается дискриминант который больше нуля или меньше нуля (колебательный или апериодический преходный процесс).

Цитата:

Mathcad 13 (символьный движок под Maple) ищет решения данного уравнения корректно (см. картинку)

У меня Mathcad 15 те же корни через solve:
413394.820656943429
-982797.21044857346926
63655.263354120156558
- видно что в 6-9 знаках после запятой они не совпадают с результатом Mathcad 13, и подстановка этих корней в уравнение не даёт точный ноль а только приблизительно (наверное достаточно чтобы догадаться что это ноль)

[GAA]

Проверю в 15 на следующей неделе. В 13, даже используя числовой (с плавающей точкой) движок, нет необходимости в замене переменной, достаточно сменить метод (LaGuerre — выбирается по умолчанию), см. прикрепление.

[DEP]

GAA в сообщении #1240317 писал(а):

Проверю в 15 на следующей неделе. В 13, даже используя числовой (с плавающей точкой) движок, нет необходимости в замене переменной, достаточно сменить метод (LaGuerre — выбирается по умолчанию), см. прикрепление.

Мне не принципиально, как говорят "инженерная точность плюс-минус 25%". Мне кажется товарищи разработчики нас не уведомляют на каком шаге итерации какого-то там численного метода поиска корня уравнения они решают завершить поиск оного корня с погрешностью не хуже или не лучше какой-то там точности представления числа с плавающей запятой. Главное начинающему, поставившему Mathcad 15 понять с кем он связался и где ноль с точностью до есть или нету.

[GAA]

svv, в некоторых случаях замена позволяет преодолеть трудности в Mathcad, но в данном случае она не приводит к улучшению результата в 13.
В прикреплении видно, что, как и ожидается, в случае символьного движка ничего не меняется, а в случае движка с плавающей запятой возня в окрестности последней значащей цифры (порядка $1\cdot 10^{-16}$).
Где-то после замены стало лучше, где-то хуже).

[svv]

Так это уже «вина» конкретных алгоритмов. Замена своё дело сделала, но если сам алгоритм неточный... (вернее, его реализация)

[GAA]

DEP, для алгебраических уравнений до четвертой степени включительно функция solve должна возвращать точные корни. В некоторых случаях эта функция может найти точные корни и для уравнений более высоких степеней. (Если коэффициенты — числа с плавающей точкой, то возвращается результат с плавающей точкой.)
Казалось бы, имея возможность получить корни с заданной точностью, можно найти приближенные значения невязки с заданной точностью, как например в Мaple (использована версия 7.0, но версия тут не важна).

> P:= z^5 -1: Z:= evalf[32](solve(P));
Z = 
1., .30901699437494742410229341718282+.95105651629515357211643933337938*I,
 -.80901699437494742410229341718282+.58778525229247312916870595463905*I, 
-.80901699437494742410229341718282-.58778525229247312916870595463905*I,
.30901699437494742410229341718282-.95105651629515357211643933337938*I
> seq(eval[32](P, z=Z[k]), k=1..5);
0., -.1e-31-.77781757615511713920653428365785e-32*I,
-.6e-31+.94876071991504312973751185717278e-31*I,
-.6e-31-.94876071991504312973751185717278e-31*I,
-.1e-31+.77781757615511713920653428365785e-32*I

Для случая коэффициентов с плавающей точкой

> Digits:= 32:
> P:= z^5 -1.0: Z:= solve(P);
Z = 1., -.80901699437494742410229341718282-.58778525229247312916870595463907*I,
-.80901699437494742410229341718282+.58778525229247312916870595463907*I,
.30901699437494742410229341718282-.95105651629515357211643933337938*I,
.30901699437494742410229341718282+.95105651629515357211643933337938*I
> seq(eval[33](P, z=Z[k]), k=1..5);
0., 0.-.13972252156122383769449788012377e-31*I,
0.+.13972252156122383769449788012377e-31*I,
-.1e-31+.77781757615511713920653428365785e-32*I,
-.1e-31-.77781757615511713920653428365785e-32*I

Однако в Mathcad 13 неожиданно для меня невязка при символьном вычислении оказывается большой. См. прикрепление.

[DEP]

GAA в сообщении #1240337 писал(а):

DEP, для алгебраических уравнений до четвертой степени включительно функция solve должна возвращать точные корни.

Я почитал справку, точность "Solve Block" контролирует системная константа CTOL, а сходимость для численных методов TOL - от неё зависит точность функции polyroots. По умолчанию они равны 0.001, но можно уменьшить.

[Евгений Машеров]

Наивный вопрос - отчего Кардано?
Если задача прикладная, и важно не освоить метод, а получить результат - как-то сразу напрашивается заведомо существующий действительный корень найти ну хоть Ньютоном, разделить на $x-x_1$ и полученное квадратное уравнение уже решать.

[GAA]

DEP, CTOL и TOL не учитываются при символьном решении уравнений (например инструмент solve на панели «Символьные» (в ru-версии)). Значения этих переменных учитываются при решении уравнений с помощью движка на базе чисел с плавающей точкой. Но легко проверить, что в рассматриваемом примере уменьшение TOL не приводит к изменению возвращаемых polyroots значений ни в 13, ни в 15.0. Однако в моей подверсии 13-ой при значительном уменьшении TOL генерируются ошибки: “This application has requested the Runtime to terminate it in unusual way”, “An internal error has occured”,... В 15.0 можно было уменьшать до 1e-15 (без видимых изменение результата).

В Mathcad 15.0 для указанных значений коэффициентов всё достаточно корректно работает. Даже глюк 13 версии, указанный в моём предыдущем сообщении, отсутствует; 15.0 выдаёт:

R(Z_1) = -4.3297218522707480238514066622814e-33+1.3965152067677095887951732637128e-32i
R(Z_2) = -8.7383899051427208767379936704027e-33-7.7871093467962764880392139402903e-33i