Форум программистов [Powered by Invision Power Board]

Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > C++ Builder > Откуда погрешность решения кубического уравнения

Автор: Нитонисе 9.1.2010, 18:25

Решаю кубическое уравнение, используя тригонометрическую формулу Виета.
Однако почему-то возникает довольно ощутимая погрешность. Это следствие несовершенства метода или весь вопрос в численных преобразованиях?
Например, при вводе коэффициентов a=b=c=d=1, естественно корень уравнения должен быть -1, однако так не получается. Почему?

Кроме того при некоторых значениях коэффициентов возникает ошибка неправильного аргумента.
Например, при a=3, b=4, c=4, d=-6, программа не может решить уравнение, потому что ей приходится вычислять натуральный логарифм, где аргументом является отрицательное число, чего быть не должно.

Код


float t = r/sqrt(pow(fabs(q),3));
float Arch = log(t+sqrt(t*t-1));
t = Arch/3;
float sh = (pow(2.7183,t)-pow(2.7183,-t))/2;
x1=-2*sqrt(fabs(q))*sh -A/3;

Однако надо отметить, что я прибегаю к вычислению натурального логарифма по той причине, что не могу сразу вычислить значения гиперболических функций Arch, Arsh, ch, sh. Эти функции не реализованы в math.h. Как можно их подключить к программе чтобы не связываться с логарифмами?

А вот сама программа.

Автор: Alexeis 9.1.2010, 19:12

Цитата(Нитонисе @ 9.1.2010, 17:25

)

довольно ощутимая погрешность

Это сколько? float гарантировано врет в 6й значащей цифре. Если значение составляет сотни, то это уже погрешность в 4м знаке после запятой.

Автор: Нитонисе 9.1.2010, 19:35

Цитата(Alexeis @ 9.1.2010, 19:12

)

Я в программе реализовал проверку погрешности.

http://www.radikal.ru

В уравнение вместо х подставляется найденный корень. Вот при решении уравнения на скриншоте - ошибка в 0.03. Я считаю это большая ошибка.

Автор: Фантом 9.1.2010, 20:15

Цитата(Нитонисе @ 9.1.2010, 19:35

)

Вот при решении уравнения на скриншоте - ошибка в 0.03. Я считаю это большая ошибка.

Не такая уж большая. На текст программы не посмотреть, но, судя по имеющемуся кусочку, никаких попыток предотвратить потерю точности Вы не предпринимали. Попробуйте заменить везде float на double и посмотреть, что получится.

Автор: Нитонисе 9.1.2010, 20:26

Вот весь текст.

Код



void __fastcall TForm1::Button1Click(TObject *Sender)
{
try
{
float a,b,c,d;
if (LEA->Text == "") a = 0;
else a = StrToFloat(LEA->Text);
if (LEA->Text == "") b = 0;
else b = StrToFloat(LEB->Text);
if (LEC->Text == "") c = 0;
else c = StrToFloat(LEC->Text);
if (LED->Text == "") d = 0;
else d = StrToFloat(LED->Text);
float A,B,C;
A = b/a;
B = c/a;
C = d/a;
float q,r,r2,q3,s;
float PI = 3.141592653589793;
q = (A*A - 3*B)/9;
r = (2*A*A*A - 9*A*B + 27*C)/54;
r2 = r*r;
q3 = q*q*q;
s = q3 - r2;
float x1,x2,x3;
int countx;
if(s > 0)
{
  float t=(acos(r/sqrt(q3)))/3;
  x1=-2*sqrt(q)*cos(t)-A/3;
  x2=-2*sqrt(q)*cos(t+2*PI/3)-A/3;
  x3=-2*sqrt(q)*cos(t-2*PI/3)-A/3;
  countx = 3;
}
if (s < 0)
{
  if (q > 0)
    {
      float t = fabs(r)/sqrt(pow(q,3));
      float Arch = log(t+sqrt(t*t-1));
      t = Arch/3;
      float znakr;
      if (r<0) znakr = -1;
      else znakr = 1;
      float ch = (pow(2.7183,t)+pow(2.7183,-t))/2;
      x1=-2*znakr*sqrt(q)*ch -A/3;
      countx = 1;
    }
  if (q < 0)
    {
      float t = r/sqrt(pow(fabs(q),3));
      float Arch = log(t+sqrt(t*t-1));
      t = Arch/3;
      float sh = (pow(2.7183,t)-pow(2.7183,-t))/2;
      x1=-2*sqrt(fabs(q))*sh -A/3;
      countx = 1;
    }
}
if (s == 0)
{
  x1 = -2*pow(r,0.3333333333333333333333) - A/3;
  x2 = x3 = pow(r,0.33333333333333333333333) - A/3;
  countx = 2;
}
if (countx == 1)
  {
    Label4->Caption = "x = " + FloatToStrF(x1,ffGeneral,5,5);
    Label5->Caption = FloatToStrF(a*x1*x1*x1+b*x1*x1+c*x1+d,ffFixed,7,5);

  }
if (countx == 2)
  {
    Label4->Caption = "x1 = " + FloatToStrF(x1,ffGeneral,5,5)
                  + "\nx2 = " + FloatToStrF(x2,ffGeneral,5,5);
    Label5->Caption = FloatToStrF(a*x1*x1*x1+b*x1*x1+c*x1+d,ffFixed,7,5)
             + "\n" + FloatToStrF(a*x2*x2*x2+b*x2*x2+c*x2+d,ffFixed,7,5);
  }
if (countx == 3)
  {
    Label4->Caption = "x1 = " + FloatToStrF(x1,ffGeneral,5,5)
                  + "\nx2 = " + FloatToStrF(x2,ffGeneral,5,5)
                  + "\nx3 = " + FloatToStrF(x3,ffGeneral,5,5);
    Label5->Caption = FloatToStrF(a*x1*x1*x1+b*x1*x1+c*x1+d,ffFixed,7,5)
             + "\n" + FloatToStrF(a*x2*x2*x2+b*x2*x2+c*x2+d,ffFixed,7,5)
             + "\n" + FloatToStrF(a*x3*x3*x3+b*x3*x3+c*x3+d,ffFixed,7,5);
  }
}
catch (...)
{
 ShowMessage("Возникли непредвиденные ошибки при работе программы.");
}

На double менял - результат аналогичен.

Автор: kemiisto 9.1.2010, 21:08

Цитата(Нитонисе @ 9.1.2010, 21:26

)

0.3333333333333333333333

1.0/3 или что-то подобное религия написать не позволяет? smile

А вообще, http://docs.sun.com/source/806-3568/ncg_goldberg.html Вам в руки. smile

Автор: Нитонисе 9.1.2010, 21:47

Цитата(kemiisto @ 9.1.2010, 21:08

)

1.0/3 или что-то подобное религия написать не позволяет?

Я писал 1/3 - это не работало, потому вставил 0.33333333333333333 smile

Добавлено через 6 минут и 12 секунд
кстати, проверял формулу Виета вручную, округляя до 20 знаков после запятой - ответ получился точно таким же.

И еще проблема с отрицательным логарифмом... как все таки можно подключить вычисление гиперболических функций?

Автор: baldina 9.1.2010, 22:13

Цитата

Код


pow(2.7183,t)

это у Вас основание натурального логарифма задано с такой поразительной точностью? И какой точности после этого Вы ожидаете в результате?

Цитата

Код


float t = r/sqrt(pow(fabs(q),3));
float Arch = log(t+sqrt(t*t-1));

Цитата

при a=3, b=4, c=4, d=-6, программа не может решить уравнение, потому что ей приходится вычислять натуральный логарифм, где аргументом является отрицательное число

при таких данных q = 3^2 - 3*4 < 0 и t должно быть fabs ® / sqrt (...)

Вобщем поглядите свой "алгоритм" еще раз внимательно, там ошибки.

Добавлено через 2 минуты и 59 секунд

Цитата

Я писал 1/3 - это не работало

1/3 и 1./3 - две большие разницы. первое - целочисленное деление, результат округляется до целого (т.е. 0)

Цитата

проверял формулу Виета вручную, округляя

чего её проверять, её 400 лет назад проверили, ошибок нет smile

а вот проверять "до 20 знака", беря e=2.7183 - сами понимаете...

Автор: Нитонисе 9.1.2010, 22:23

Цитата(baldina @ 9.1.2010, 22:13

)

Указанный фрагмент относится к вычисление гиперболического косинуса, который равен:
ch(x) = [e^x + e^(-x)]/2
Я взял е равным 2.7183. А чему равно е?

Цитата(baldina @ 9.1.2010, 22:13

)

при таких данных q = 3^2 - 3*4 < 0 и t должно быть fabs r / sqrt (...)

http://ru.wikipedia.org/wiki/Тригонометрическая_формула_Виета Посмотрите ссылку, там все же r берется не по модулю.

Цитата(baldina @ 9.1.2010, 22:13

)

Вобщем поглядите свой "алгоритм" еще раз внимательно, там ошибки.

А где еще?

Автор: Фантом 9.1.2010, 22:28

Цитата(Нитонисе @ 9.1.2010, 22:23

)

Я взял е равным 2.7183. А чему равно е?

А зачем вообще его брать? Функция exp в math.h есть.

Автор: Нитонисе 9.1.2010, 22:43

Цитата(Фантом @ 9.1.2010, 22:28

)

А зачем вообще его брать? Функция exp в math.h есть.

Вот так объявлена функция exp
double exp (double __x)
И что мне передать в нее в качестве параметра х?

Автор: Alexeis 9.1.2010, 22:48

Нитонисе
M_E -- The base of natural logarithms (e).

http://msdn.microsoft.com/en-us/library/4hwaceh6.aspx

Автор: Фантом 9.1.2010, 22:54

Цитата(Нитонисе @ 9.1.2010, 22:43

)

И что мне передать в нее в качестве параметра х?

Показатель экспоненты, надо думать. smile

Автор: Нитонисе 9.1.2010, 23:16

Цитата(Alexeis @ 9.1.2010, 22:48

)

M_E -- The base of natural logarithms (e).

Заменил экспоненту на M_E, число пи на M_PI. Результат тот же.

Добавлено через 1 минуту и 26 секунд
Так а что насчет гиперболических функций? Как их подключить?

Добавлено через 2 минуты и 4 секунды

Цитата(Фантом @ 9.1.2010, 22:54

)

Показатель экспоненты, надо думать.

Ну и какой у нас показатель экспоненты?