Решение кубического уравнения
Пост для тех, кто, как и я, не понял этого в школе ☺
Из школы осталось воспоминание, что квадратное уравнение решить проще простого - дискриминант там, больше нуля, значит, есть действительные корни, и т.д. Да, кстати, о комплексных корнях и о кратных говорить не будем, считаем, что у кубического уравнения есть три действительных корня, и их нужно найти. Бывает, что уравнение четвёртого порядка удаётся свести к биквадратному, там тоже просто. А как быть с кубическим?
Здесь школьное воспоминание говорило о том, что все три корня придётся искать численно. А это значит, что надо найти подбором три отрезка, на которых левая часть уравнения меняет знак, и в каждом из отрезков запустить поиск решения, например, половинными делениями. На самом деле не всё так уныло. Рассмотрим уравнение:
ax3+bx2+cx+d=0
Сразу решить мы его не можем. Но мы довольно легко можем определить локальные максимум и минимум на этой кривой. Казалось бы, они нам совершенно не нужны. Но, если корня три, то один из них окажется где-то между этими локальными экстремумами. Поэтому считаем производную, приравниваем её нулю и решаем получившееся квадратное уравнение:
3ax2+2bx+c=0
Находим его корни, обозначим их x'1 и x'2. Где-то между ними лежит один из корней уравнения. Ищем этот корень половинными делениями. В приведённом ниже коде вместо иксов со штрихами написаны лямбды, но смысла это не меняет.
import math
a = 1
b = -46.50489
c = 400.4657
d = -483.3768
def myfunc(lamb):
return a * lamb ** 3 + b * lamb ** 2 + c * lamb + d
eps = 1E-12
lamb1 = 5.1666
lamb2 = 25.8366
lamb3 = (lamb1+lamb2)/2
while abs(myfunc(lamb3)) > eps:
f1 = myfunc(lamb1)
f2 = myfunc(lamb2)
f3 = myfunc(lamb3)
if math.copysign(1.0, f1) == math.copysign(1.0, f3):
lamb1 = lamb3
else:
lamb2 = lamb3
lamb3 = (lamb1+lamb2)/2
print(lamb3) Теперь у нас есть один корень уравнения, назовём его x1. А это значит, что от исходного кубического уравнения мы можем перейти к квадратному, вынося (x-x1) за скобки, и решить его обычным способом:
(x-x1)(ax2+(ax1+b)x+(ax12+bx1+c))=0
Из школы осталось воспоминание, что квадратное уравнение решить проще простого - дискриминант там, больше нуля, значит, есть действительные корни, и т.д. Да, кстати, о комплексных корнях и о кратных говорить не будем, считаем, что у кубического уравнения есть три действительных корня, и их нужно найти. Бывает, что уравнение четвёртого порядка удаётся свести к биквадратному, там тоже просто. А как быть с кубическим?
Здесь школьное воспоминание говорило о том, что все три корня придётся искать численно. А это значит, что надо найти подбором три отрезка, на которых левая часть уравнения меняет знак, и в каждом из отрезков запустить поиск решения, например, половинными делениями. На самом деле не всё так уныло. Рассмотрим уравнение:
ax3+bx2+cx+d=0
Сразу решить мы его не можем. Но мы довольно легко можем определить локальные максимум и минимум на этой кривой. Казалось бы, они нам совершенно не нужны. Но, если корня три, то один из них окажется где-то между этими локальными экстремумами. Поэтому считаем производную, приравниваем её нулю и решаем получившееся квадратное уравнение:
3ax2+2bx+c=0
Находим его корни, обозначим их x'1 и x'2. Где-то между ними лежит один из корней уравнения. Ищем этот корень половинными делениями. В приведённом ниже коде вместо иксов со штрихами написаны лямбды, но смысла это не меняет.
import math
a = 1
b = -46.50489
c = 400.4657
d = -483.3768
def myfunc(lamb):
return a * lamb ** 3 + b * lamb ** 2 + c * lamb + d
eps = 1E-12
lamb1 = 5.1666
lamb2 = 25.8366
lamb3 = (lamb1+lamb2)/2
while abs(myfunc(lamb3)) > eps:
f1 = myfunc(lamb1)
f2 = myfunc(lamb2)
f3 = myfunc(lamb3)
if math.copysign(1.0, f1) == math.copysign(1.0, f3):
lamb1 = lamb3
else:
lamb2 = lamb3
lamb3 = (lamb1+lamb2)/2
print(lamb3) Теперь у нас есть один корень уравнения, назовём его x1. А это значит, что от исходного кубического уравнения мы можем перейти к квадратному, вынося (x-x1) за скобки, и решить его обычным способом:
(x-x1)(ax2+(ax1+b)x+(ax12+bx1+c))=0