Решение кубического уравнения (HP-35s)

Submitted by basvic on Tue, 28/10/2014 - 12:34

Решение невырожденного кубического уравнения с вещественными коэффициентами методом Кардано. Поскольку этот калькулятор не умеет решать такие уравнения из коробки, я посочинял немного. Возможно пригодится кому-либо.

Входные данные - уравнение вида ax3+bx2+cx+d=0

Коэффициенты вводятся по запросу калькулятора c нажатием клавиши R/S

Выходные данные - корни (вещественные или комплексные) находятся в трёх нижних регистрах стека x3->Z, x2->Y, x1->X

Я тут окинул взором, что получилось и решил сократить код на 23 команды выделив 2 подпрограммы...

C001	LBL C
C002	RPN
C003	RAD
C004	SF 10
C005	AX^3+BX^2+CX+D=0
C006	PSE
C007	INPUT A 
C008	INPUT B
C009	INPUT C
C010	INPUT D

C011	RCL B       //поправка B/(3A)
C012	RCL A
C013	/
C014	3
C015	/
C016	STO E

C017	3           //Вычисление P
C018	RCL C
C019	*
C020	RCL A
C021	*
C022	RCL B
C023	x2
C024	-
C025	RCL A
C026	x2
C027	3
C028	*
C029	/
C030	STO P

C031	RCL D       //Вычисление Q
C032	RCL A
C033	/
C034	RCL B
C035	RCL C
C036	*
C037	3
C038	/
C039	RCL A
C040	x2
C041	/
C042	-
C043	RCL B
C044	3
C045	yX
C046	2
C047	*
C048	27
C049	/
C050	RCL A
C051	3
C052	yX
C053	/
C054	+
C055	STO Q

C056	x2         //Дискриминант
C057	4
C058	/
C059	RCL P
C060	3
C061	yX
C062	27
C063	/
C064	+
C065	STO S

C066	x>0?          //Выбор ветви решения
C067	GTO C096
C068	x<0?
C069	GTO C148

C070	RCL Q         // Дискриминант = 0 Корни вещественные кратные
C071	+/-
C072	2
C073	/
C074	3
C075	x SQR y
C076	STO U
C077	RCL B
C078	RCL A
C079	/
C080	3
C081	/
C082	STO V
C083	2
C084	RCL U
C085	*
C086	RCL V
C087	-
C088	STO X
C089	RCL U
C090	+/-
C091	RCL V
C092	-
C093	STO Y
C094	STO Z
C095	GTO C193

C096	RCL S           // Дискриминант > 0. Один вещественный корень и два комплексных
C097	SQRT x
C098	STO T
C099	RCL Q
C100	2
C101	/
C102	+/-
C103	STO W
C104	RCL T
C105	+
C106	3
C107	x SQR y
C108	STO U
C109	RCL W
C110	RCL T
C111	-
C112	3
C113	x SQR y
C114	STO V
C115	RCL U
C116	+
C117	STO R
C118	RCL U
C119	RCL V
C120	-
C121	STO I
C122	RCL R
C123	RCL E
C124	-
C125	STO X
C126	XEQ C133
C127	-
C128	STO Y
C129	XEQ C133
C130	+
C131	STO Z
C132	GTO C193

C133	RCL R         // подпрограмма 1
C134	2
C135	/
C136	+/-
C137	RCL E
C138	-
C139	RCL I
C140	3
C141	SQRT x
C142	*
C143	2
C144	/
C145	i
C146	*
C147	RTN

C148	RCL Q       // Дискриминант < 0 - три разных вещественных корня
C149	2
C150	/
C151	+/-
C152	RCL P
C153	+/-
C154	3
C155	yX
C156	27
C157	/
C158	SQRT x
C159	/
C160	ACOS
C161	STO F
C162	2
C163	RCL P
C164	+/-
C165	3
C166	/
C167	SQRT x
C168	*
C169	STO U
C170	0
C171	XEQ C180
C172	STO X
C173	2
C174	XEQ C180
C175	STO Y
C176	-2
C177	XEQ C180
C178	STO Z
C179	GTO C193

C180	RCL F        // Подпрограмма 2
C181	x<->y
C182	pi
C183	*
C184	+
C185	3
C186	/
C187	COS
C188	RCL U
C189	*
C190	RCL E
C191	-
C192	RTN

C193	VIEW X       // Вывод результатов
C194	VIEW Y
C195	VIEW Z
C196	STOP
C197	GTO C001

Символы отличные от изображённых на клавиатуре калькулятора
* - умножение
/ - деление
SQRT x - квадратный корень
x SQR y - корень степени X из содержимого Y регистров стека
x<->y - обмен данными между регистрами стека X и Y
pi - 3,1415927....

Метки публикаций: 
Кубическое уравнение
Алгебра
HP (Hewllet Packard)

Comments

st's picture
Submitted by st on Wed, 29/10/2014 - 11:17

Мощно, от листинга повеяло старой школой :) Наверняка программу можно подсократить. В "НиЖ" лет 30 назад обсуждался этот вопрос: надо ли "полировать код" или "считает - и хорошо".

Submitted by basvic on Wed, 29/10/2014 - 16:10

Сейчас мне очевидно как сократить ещё на 5 команд, но это уже будет во вред удобочитаемости кода (если об этом языке можно так говорить). Поэтому я остановился на таком варианте. Ну а чтобы выделить подпрограммы, сначала надо было отладить код в простейшем варианте, а затем взглянуть на весь код целиком. Чем короче код, тем меньше шансов совершить ошибку при вводе.
Ну а поскольку программа рассчитана на то, что ею будут пользоваться другие люди, то смысл полировать код есть.

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Submitted by Vlad_AB on Wed, 29/10/2014 - 17:18

Для решения кубического уравнения можно использовать программку: 

F001 LBL F
F002 A*X^3+B*X^2+C*X+D
F003 RTN

(LN=26)

E001 LBL E
E002 RPN
E003 SF 10
E004 AX^3+Bx^2+CX+D=0
E005 PSE
E006 INPUT A
E007 INPUT B
E008 INPUT C
E009 INPUT D
E010 CF 10
E011 FN= F
E012 SOLVE X
E013 GTO E015
E014 RTN
E015 A->E
E016 B+A*X->F
E017 C+B*X+A*X^2->G
E018 (-F-((SQ(F)-4*E*G)+i*0)/(2*E)->Y
E019 (-F+((SQ(F)-4*E*G)+i*0)/(2*E)->Z
E020 VIEW X
E021 VIEW Y
E022 VIEW Z
E023 RTN

(LN=180)

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

Submitted by basvic on Wed, 29/10/2014 - 22:13

Обязательно попробую Ваше альтернативное решение
Кстати, обратил внимание, что в строках 18 и 19 на 5 открытий скобок только 4 закрытия

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Submitted by basvic on Wed, 29/10/2014 - 23:24

Потестировал я Ваш альтернативный вариант дописав скобку в конец выражения и вот что получилось:
x3+3x2-6x-8=0 получилось (-4 -3.5+i0 5.5+i0), а должно быть (-1 2 -4)
2x3-11x2+12x+9=0 получилось (-0.5 12+i0 12+i0), а должно быть (-0,5 3 3)
x3+x2+6x+16=0 получилось (-2 16.5+i0 -14.5+i0), а должно быть (-2 0.5+i*2.7839... 0.5-i*2.7839...)
Один вещественный корень находит, но два других - увы. Пробовал и с другим положением упущенной скобки но результат тоже далёк от истины. Ваша версия нуждается в отладке.

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Submitted by Vlad_AB on Thu, 30/10/2014 - 04:59

Потестировал я Ваш альтернативный вариант дописав скобку в конец выражения и вот что получилось:
x3+3x2-6x-8=0 получилось (-4 -3.5+i0 5.5+i0), а должно быть (-1 2 -4)
2x3-11x2+12x+9=0 получилось (-0.5 12+i0 12+i0), а должно быть (-0,5 3 3)
x3+x2+6x+16=0 получилось (-2 16.5+i0 -14.5+i0), а должно быть (-2 0.5+i*2.7839... 0.5-i*2.7839...)
Один вещественный корень находит, но два других - увы. Пробовал и с другим положением упущенной скобки но результат тоже далёк от истины. Ваша версия нуждается в отладке.

:-) Отладка была проведена перед записью комментария с программкой. В моём кальке всё считает верно. К сожалению сразу не посмотрел правильно ли набран текст этой программки в опубликованном комментарии. Увы, после записи вашего комментария редактирование моего комментария становится недоступным.

Текст программки: 

F001 LBL F
F002 A*X^3+B*X^2+C*X+D
F003 RTN

(LN=26)

E001 LBL E
E002 RPN
E003 SF 10
E004 AX^3+Bx^2+CX+D=0
E005 PSE
E006 INPUT A
E007 INPUT B
E008 INPUT C
E009 INPUT D
E010 CF 10
E011 FN= F
E012 SOLVE X
E013 GTO E015
E014 RTN
E015 A->E
E016 B+A*X->F
E017 C+B*X+A*X^2->G
E018 (-F-((SQ(F)-4*E*G)+i*0)^0.5)/(2*E)->Y
E019 (-F+((SQ(F)-4*E*G)+i*0)^0.5)/(2*E)->Z
E020 VIEW X
E021 VIEW Y
E022 VIEW Z
E023 RTN

(LN=180)

Очевидные недостатки программы:

  1. Скорость работы.
  2. Для некоторых уравнений SOLVE не может найти корень.
  3. Для последних кратных действительных корней приписывает мнимую единицу с нулём.

Скорость никак не изменишь, приписывание мнимой единицы является особенностью работы калькулятора с комплексными числами, а вот поиск корня функцией SOLVE можно улучшить предварительным заданием начальных приближений.

Submitted by basvic on Thu, 30/10/2014 - 10:50

Теперь работает. Однако использование численных методов в случае, когда есть аналитическое решение - это с точки зрения математики моветон. Хотя, конечно, когда известен один корень, то, можно поделить многочлен 3-й степени на многочлен первой и получить квадратное уравнение.

Это напоминает историю с зеркалом Хаббла. Одно (запасное) сделали классическим методом, второе альтернативным и вывели на орбиту зеркало с ошибкой фигуризации поверхности.

Как быть, если в ответственный момент программа не найдёт численным методом первый корень? Не много людей держат в голове метод Кардано, чтобы решить уравнение ручным счётом.

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Submitted by Vlad_AB on Tue, 04/11/2014 - 18:03

Теперь работает. Однако использование численных методов в случае, когда есть аналитическое решение - это с точки зрения математики моветон. Хотя, конечно, когда известен один корень, то, можно поделить многочлен 3-й степени на многочлен первой и получить квадратное уравнение.

Моветон – это использование в программе только чистого RPN, когда можно использовать и уравнения, особенно когда речь идёт о читабельности кода.

Это напоминает историю с зеркалом Хаббла. Одно (запасное) сделали классическим методом, второе альтернативным и вывели на орбиту зеркало с ошибкой фигуризации поверхности.

Пример с телескопом Хаббла интересный, но некорректный в той же степени, как и «вывод» о японской бензопиле (см. анекдот далее). В обоих случаях технология не подводила, а свою негативную роль сыграл человеческий фактор. Если я правильно помню, то ошибка в обработке поверхности зеркала возникла из-за шайбы, а не из-за альтернативности технологии (бензопила сломалась из-за…).

Лесорубам подарили новую японскую бензопилу. Подставили доску: - Вжик! - сказала японская бензопила. - Ух ты! - сказали лесорубы. Подставили бревно: - Вжик! - сказала японская бензопила. - Ух ты! - сказали лесорубы. Подставили железный лом: - Крррр....! - сказала японская бензопила. - Аaaaaa, бл*! - сказали лесорубы.

Как быть, если в ответственный момент программа не найдёт численным методом первый корень? Не много людей держат в голове метод Кардано, чтобы решить уравнение ручным счётом.

Функция, определяемая кубическим уравнением, если судить по руководству пользователя является примером хорошей функции, для которой функция SOLVE может найти корень. А в случае какого-либо «сбоя» программу всегда можно дополнить автоматическим перезаданием начального приближения для корня (об этом писал в предыдущем комментарии). Действительно, вашу программу с формулой Кардано сложнее запомнить и набрать, чем альтернативную программку.

Submitted by basvic on Wed, 05/11/2014 - 12:12

Ай-ай... Пошёл посыпать голову пеплом...

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Submitted by Vlad_AB on Sun, 09/11/2014 - 07:36

Написал программку с аналитическим алгоритмом решения кубического уравнения с использованием уравнений (по формуле Кардано):

K001 LBL K
K002 CF 10
K003 INPUT A
K004 INPUT B
K005 INPUT C
K006 INPUT D
K007 B/3/A->E
K008 C/A-3*SQ(E)->P
K009 2*E^3-C*E/A+D/A->Q
K010 (P^3/27+Q^2/4)^0.5i0->T
K011 (-Q/2+T)^(1/3)->R
K012 ABS
K013 x=0?
K014 (-Q/2-T)^(1/3)->R
K015 -P/3/R->S
K016 R+S-E->X
K017 ABS(X)*COS(ARG(X))->X
K018 A->E
K019 B+A*X->F
K020 C+F*X->G
K021 -F/2/E->I
K022 (I^2-G/E)^0.5i0->J
K023 -
K024 STO Y
K025 I+J->Z
K026 VIEW X
K027 VIEW Y
K028 VIEW Z
K029 CLSTK
K030 RTN

(LN=263)

На всех возможных подслучаях не тестировал.