Помощь - Поиск - Пользователи - Календарь
Полная версия: Вычисление относительной погрешности для функции.
Форум «Всё о Паскале» > Delphi, Assembler и другие языки. > Другие языки
Krjuger
В общем это некое продолжение моей прошлой темы,только обрастающее новыми подробностями.
У меня дана функция x(expx-1).Я эту функцию раскладываю я ряд Тейлора и получаю сумма от 1 до N от xn+1/n!.
Чтобы найти относительную погрешность мне надо а -n ый член разделить на сумму ряда,все это по модулю.Ну в общем то ,я думаю, вы и так понимаете как это делать.Суть заключается в том,у меня опять есть ограниченная разрядность мантисы и мне надо посмотреть как она будет влиять на результат.И как будет меняться N при которых мы будем выходить за граници возможностей нашей машины.

Я немного абстрагировался от этой задачи.что я сделал, при помощи маткада я посчитал,при каком N будет достигаться относительная погрешность на грани машинного эпсилон 10-16,это N =24.Да забыл сказать,что точка в которой мы раскладываем ряд у меня дана.Это -2.3.Затем я высчитал значение этой погрешности она составила 8.602898672363349*10^-16.Дальше я посомтрел как оно себя будет вести при разрядах мантисы от 10 до 25,при фиксированных исходных данных.Для рязрядности с 10 до 16 я получил ожидаемый результат,но дальше начало твориться что то вообще непонятное.Ну или точнее я не могу понять, как это интерпретировать.
ИзображениеФайл с самой программой тоже прикрепляю.
Krjuger
Я не много изменил основную программу,она стала выглядить так

void main()
{
double num=0;
double res=0;
double x1,x2;
int m=0;
int N;
x1=-2.3;
x2=-17;
for (int i=5;i<=20;i++)
{
cout << "Razryadnost* manticy:" << i << endl;
N=0;
for (int s=1;s<25;s++)
{
res=DeltaPogr(x1,s,i);
if (res!=0)
{
cout << setprecision (10) << res << endl;
N++;
}
}
cout << "Kol-vo chlenov ryada:" << N << endl;
}
}


Чтобы стало немного нагляднее.Принципи,большинство цифр получились вполне ожидаемыми,но возникли вопросы,почему,при разрядности мантисы 19 получается самое точное число ,а именно 8.603*10^-16,но дальше начинается уже какой то бред,так же непонятно,потому при 19,по идеи это должно было произойти при 16,либо я что то сильно недопонимаю.И я не понял почему при 20 я не получил на 1 значащую цифру более точный результат?В общем вопросов много....
volvo
Что-то тема заскучала... Давай сделаем вброс небольшой.

1)
Цитата
что я сделал, при помощи маткада я посчитал,при каком N будет достигаться относительная погрешность на грани машинного эпсилон 10-16,это N =24.
Неверно. Надо проверять не на грани машинного эпсилон, а на грани эпсилон для твоего типа, мантисса которого содержит ограниченное число разрядов. Для m = 5 и m = 11 это будут разные значения, и нет смысла вычислять с точностью 10-16, если m = 5. Максимальная точность - только для максимальных длин мантиссы, и то не в типе float, а в long double. Соответственно, вычислять надо не N первых членов ряда, а пока очередной член не будет меньше Epsilon для типа.

2)
Цитата
Дальше я посомтрел как оно себя будет вести при разрядах мантисы от 10 до 25,при фиксированных исходных данных.
На всех значениях от 10 до 25 посмотреть не могу, а вот для M = 5, M = 8 и M = 11 - проверил. Извиняй, в Сях нет соответствующих средств, поэтому проверил там, где есть:

with Ada.Text_IO; use Ada.Text_IO;

procedure Main is

generic
type Float_Type is digits <>;
procedure Check(X : Float_Type);

procedure Check(X : Float_Type) is
s, nom, nxt, last : Float_Type;
n, denom : Long_Long_Integer := 1;
begin
Ada.Text_IO.Put_Line ("Mantissa : " & Integer'Image (Float_Type'Mantissa));
Ada.Text_IO.Put_Line ("Epsilon : " & Float'Image (Float_Type'Epsilon));

s := 0.0;
nom := X * X;
loop
nxt := nom / Float_Type (denom);
Ada.Text_IO.Put_Line ("nxt = " & Float_Type'Image(nxt) &
" n! = " & Long_Long_Integer'Image(denom));
exit when abs(nxt) < Float_Type'Epsilon;
s := s + nxt;
nom := nom * X;
n := n + 1;
denom := denom * n;
last := nxt;
end loop;
Ada.Text_IO.Put_Line ("Last item : " & Float_Type'Image(last));
Ada.Text_IO.Put_Line ("S : " & Float_Type'Image(s));
Ada.Text_IO.Put_Line ("N : " & Long_Long_Integer'Image(n));

Ada.Text_IO.Put_Line ("Relative error : " &
Float'Image( abs(Float(last) / Float(s)) ));
Ada.Text_IO.New_Line;
end Check;

type MyFloat_05 is new Float digits 1;
procedure Check_05 is new Check(MyFloat_05);

type MyFloat_08 is new Float digits 2;
procedure Check_08 is new Check(MyFloat_08);

type MyFloat_11 is new Float digits 3;
procedure Check_11 is new Check(MyFloat_11);


begin
Check_05 (X => -2.3);
Check_08 (X => -2.3);
Check_11 (X => -2.3);
end Main;
Результат:
Mantissa :  5
Epsilon : 6.25000E-02
nxt = 5.3E+00 n! = 1
nxt = -6.1E+00 n! = 2
nxt = 4.7E+00 n! = 6
nxt = -2.7E+00 n! = 24
nxt = 1.2E+00 n! = 120
nxt = -4.7E-01 n! = 720
nxt = 1.6E-01 n! = 5040
nxt = -4.5E-02 n! = 40320
Last item : 1.6E-01
S : 2.1E+00
N : 8
Relative error : 7.38203E-02

Mantissa : 8
Epsilon : 7.81250E-03
nxt = 5.3E+00 n! = 1
nxt = -6.1E+00 n! = 2
nxt = 4.7E+00 n! = 6
nxt = -2.7E+00 n! = 24
nxt = 1.2E+00 n! = 120
nxt = -4.7E-01 n! = 720
nxt = 1.6E-01 n! = 5040
nxt = -4.5E-02 n! = 40320
nxt = 1.1E-02 n! = 362880
nxt = -2.6E-03 n! = 3628800
Last item : 1.1E-02
S : 2.1E+00
N : 10
Relative error : 5.51081E-03

Mantissa : 11
Epsilon : 9.76563E-04
nxt = 5.29E+00 n! = 1
nxt = -6.08E+00 n! = 2
nxt = 4.66E+00 n! = 6
nxt = -2.68E+00 n! = 24
nxt = 1.23E+00 n! = 120
nxt = -4.73E-01 n! = 720
nxt = 1.55E-01 n! = 5040
nxt = -4.47E-02 n! = 40320
nxt = 1.14E-02 n! = 362880
nxt = -2.63E-03 n! = 3628800
nxt = 5.49E-04 n! = 39916800
Last item : -2.63E-03
S : 2.07E+00
N : 11
Relative error : 1.26910E-03

[2011-04-15 12:00:33] process terminated successfully (elapsed time: 00.23s)

Могу проверить еще и для остальных (M = 15, 18, 21, ...), если надо.
Krjuger
Я наверно со совсем корректно написал то,что сделал.
Цитата

Неверно. Надо проверять не на грани машинного эпсилон, а на грани эпсилон для твоего типа, мантисса которого содержит ограниченное число разрядов. Для m = 5 и m = 11 это будут разные значения, и нет смысла вычислять с точностью 10-16, если m = 5. Максимальная точность - только для максимальных длин мантиссы, и то не в типе float, а в long double. Соответственно, вычислять надо не N первых членов ряда, а пока очередной член не будет меньше Epsilon для типа.

Да я это понимаю.Я их не сравнивал с машинной эпсилон.
В общем,лень мое горе,напишу полностью,как есть информацию.Изначально у меня было 2 задачи 1 нужно было сделать в маткаде,а 2 на языке.
Цитата

Задача 1.2 Используя разложение заданной функции F(x) в ряд Тейлора в окрестности нуля, вычислить значения функции в двух точках x1и x2. Вычисление частичных сумм ряда производить до тех пор, пока отношение прибавляемого члена к частичной сумме не станет меньше машинного эпсилон..
ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ:
1. Разложить заданную функцию F(x) в ряд Тейлора в окрестности нуля.
2. Выписать величину относительной погрешности, вычисляя ее как отношение прибавляемого члена к накопленной частичной сумме , взятое по модулю.
3. Определить количество членов ряда , при котором величина относительной погрешности станет
меньше машинного эпсилон.

Это я сделал сдал и бог бы с ним.В маткаде было выставлено считать до 16 цифр,результат я показал ранее.
Цитата
Задача 1.4. Составить программу на алгоритмическом языке, моделирующую вычисления на ЭВМ с ограниченной разрядностью m. Решить задачу 1.2 , используя эту программу. Составить график зависимости относительной погрешности от количества разрядов m= 4,5,…8.

Вот собственно то,что я пытаюсь сотворить.
Я немного отошел от задания,а именно сделал проверку при мантисе не 4,5,…8,а 10-25,чтобы посмотреть поведение.В моей голове уложилось такое мнение,что если мантиса имеет n чисел,то для этой мантисы машинное эпсилон будет равно 10-n,Так же я понимаю,что существует машинная эпсилон для тех стантартных типов данных,что есть я языке.В маткаде при 16 цифрах,я получал результат 8.602898672363349*10^-16,что вполне логично,это было границей точности,так как следующее число было уже 10-17,что выходит за рамки.Все это действо получалось при суммировании 24 членов ряда.Я думаю вполне логично,что при мантисе от 10 до 25 где нить должно будет появиться это число.Когда я выставляю мантису 16 и суммирую 24 раза я получаю вполне ожидаемый результат 8.602898672363349*10^-16 превращается в 9*10^-16,что вполне логично.Эта логика сохраняется до мантисы в 19.Происходить просто увеличение на 1 цифру точности числа,при 19 я получаю 8.603*10^-16,но при 20 происходит ,то чего я не ожидал.Вместо улучшения еще на 1 цифру я получаю совершенно другие результаты.Я пробовал менять double на long double,думая,что упираюсь в его точность,но меня ждал,сурприз,ни в одном тесте цифры не изменились.

Тип float я в своей программе нигде и не использовал.Собственно,вопросы которыми я задаюсь следующие.Почему фигня начинается при мантисе 19,а не после 16,я спрашивал у преподавателя,он что то невнятное сказал,про то ,что когда сопроцессор считает,он использует расширенную точность,меня это как то невоодушевило,слишком расплывчато.
В ваш код я очень старался вникнуть,но увы,с адой мое знакомство нулевое,так что я понял достаточно мало.Единственное я понял,что вы высчитывали факториал и слешующий элемент последовательности.

Я не понял почему вы вызывали для мантисы 5 type MyFloat_05 is new Float digits 1;
И в результате получали Last item : 1.6E-01,тобиш 0.16,вроде как мантиса 5 должна упереться в 10-5,либо я уже совсем ничего не понимаю.Кстати вы уж извините,но я руками посчитал несколько первых цленов ряда,и с вашими цифрами они никак не совпадают.

Так что я прошу прокоментировать ваш код,ну или сделать небольшой вброс теории, которой вы руководствовались.Ну и если я где то что то смутно написал,то спросить,чтобы я прямо ответил.


Да я забыл сказать,Для второй точки там число суммирований членов ряда при 16 получается около 70,так что я не стал реализовывать факториал рекурсивно,такая махина ни в один стандартный тип не влезет,использовать длинную арифметику тоже не стал,на мой взгляд не стоит овчинка выделки,да и как сделать это не знаю,не сталкивался,так что я поступил следующим образом,переходя в цикле на очередной член ряда я его в цикле делил на число от 2 до номера текущего члена ряда,и каждый раз при делении укорачивал подстраиваясь под мантису.Может это создаст еще некоторую погрешность,но мне легче обосновать преподавателю эту вынужденную погрешность,чем избавиться от нее.
volvo
Цитата
Я не понял почему вы вызывали для мантисы 5 type MyFloat_05 is new Float digits 1;
Я не для какой-то мантиссы делал специфическое значение digits, я просто создал три новых типа, каждый из которых ограничивает вычисления определенной точностью. Первый - одной десятичной цифрой после запятой, второй - двумя, третий, соответственно, тремя. Зачем? А затем, чтобы все эти 3 типа имели разное представление: разную длину мантиссы, в данном случае. При желании можно сделать еще и длину экспоненты разной, но меня это в данном случае меньше интересовало. Потом проверил, что для сконструированного типа возвращает атрибут 'Mantissa - это как раз и будет число бит в мантиссе (как ты на 5 битах собрался достичь точности 10-5 для меня - загадка). Атрибут 'Epsilon возвращает ту самую точность, до достижения которой надо продолжать вычисления (здесь мне ничего вручную делать не надо - в отличие от тебя, творящего непонятно что с числами - в случае Ады это забота компилятора. Если я сказал, что мне достаточно одной цифры после запятой, и 'Mantissa показала, что M=5, то все вычисления будут производиться именно так, как они и должны производиться при пятибитной мантиссе).

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

Цитата
Собственно,вопросы которыми я задаюсь следующие.Почему фигня начинается при мантисе 19,а не после 16
Собственно, вопрос: а почему фигня должна начинаться при длине мантиссы = 16? И с какой стати ты решил, что те результаты, которые ты получил - фигня? Я, например, вообще ничего не вижу на твоем скриншоте. Выведи хотя бы то же самое, что выводил я: длина мантиссы, последний член ряда, сумма ряда. Вот тогда и посмотрим, что у тебя творится...

Цитата
Кстати вы уж извините,но я руками посчитал несколько первых членов ряда,и с вашими цифрами они никак не совпадают.
Вы у ж меня тоже извините, но я посчитал сумму ряда (оно же - значение x(ex - 1)) при заданном X. Что характерно - сумма (при использовании точных типов) совпадает идеально. Как может быть, что член ряда не совпадает, а их сумма - таки да? И... Что именно не совпадает?

Первый: (-2.3)2 / 1! = 5.29 (у меня 5.3E+00)
Второй: (-2.3)3 / 2! = -6.0835 (у меня -6.1E+00)

Дальше продолжать? С учетом точности представления (не забыл? 5 бит мантиссы - это 1 десятичная цифра после запятой) как раз все совпадает...

В общем, убеждаюсь в очередной раз: тебе лучше не отвечать, ибо я ж еще и виноват в том, что ты не умеешь читать и считать. dry.gif Удачи тебе в дальнейших обвинениях других...

Вот вычисления при длине мантиссы 15, 18, 21, 25, 28 и 31 бит (Показать/Скрыть)
может быть пригодится...
Krjuger
Почему,реакция сразу негативная,я же не говорил,что вы не правы,я лиш сказал,что мои расчеты и ваши отличаются, именно поэтому я попросил привести аргументы,которыми вы руководствовались,чтобы найти ошибку.
Насчет всего выше сказанного буду пытаться обмозговать.

Скриншот исправлю,он относился к первой редакции,до изменений,дальше стано намного читабельней.

Я нашел почему получались разные результаты,да,это была моя ошибка,ну а точнее неточность,я выводил сразу относительную погрешность,а не а-n-ый член рядаи даже не сумму.Так,по числам у меня воспадает с вашими.
Насчет эпсилон,я покопался в рукописах своих,и ума не приложу,почему меня так переклинило на этой точности.Буду работать ,чуть позже выложу новый вариант.
Krjuger
В общем,я нашел основные ошибки и исправил их.Вышло то,что и получалось у тебя на скринах,volvo.
Я выложу код,который считает an,чтобы можно было сравнить с твоими результатами.

#include <iostream>
#include <math.h>
#include <iomanip>


using namespace std;

double Round(long double chislo, int n)
{
long double result;
__int64 iChislo;
double drob;
int k=0;
if (chislo<0)
k=-1;
else
k=1;
iChislo = (__int64) chislo;
drob = chislo - iChislo;
result=((__int64) (drob*pow(10.0,n)+k*0.5))/pow(10.0,n);
result = iChislo + result;
return result;
}


double Cut(long double value,int m)
{
int n=0;
long double temp=abs(value);
while (temp>=1)
{
temp=temp/10;
n++;
}
value=Round(value,abs(m-n));
return value;
}

double Sum(long double num,int n,int m)
{
long double temp=0;
long double tmp=0;
for (int i=1;i<=n;i++)
{
tmp=pow(num,i+1);
tmp=Cut(tmp,m);
for (int j=2;j<=i;j++)
{
tmp=tmp/j;
tmp=Cut(tmp,m);
}
temp+=tmp;
}
return temp;
}

double DeltaPogr (long double x,int n,int m)
{
long double temp=0;
long double an=0;
an=pow(x,n+1);
an=Cut(an,m);
for (int j=2;j<=n;j++)
{
an=an/j;
an=Cut(an,m);
}
temp=abs(an/Sum(x,n,m));
temp=Cut(temp,m);
return an;
}

void main()
{
long double num=0.0;
long double res=0.0;
long double x1,x2;
long double eps=0.0;
long double temp=0.0;
int m=0;
int s=0;
int N;
x1=-2.3;
x2=-17;
for (int i=5;i<=8;i++)
{
cout << "Razryadnost* manticy:" << i << endl;
N=0;
m=0;
s=0;
eps=1/pow(2.0,i-1);
temp=eps;
while (temp<=1)
{
temp*=10;
m++;
}
cout << "Epsilon:"<< eps << endl;
res=DeltaPogr(x1,1,m+1);
while (abs(res)>=eps)
{
s++;
res=DeltaPogr(x1,s,m+1);
cout << setprecision (10) << res << endl;
N++;
}
cout << "Kol-vo chlenov ryada:" << N << endl;
}
}


Пока что написано достаточно коряво,много повторений практически идентичных кусков кода.
Изображение
Krjuger
Я немного поработал на округлением,а то до этого не совсем корректно выводились данные и в итоге пришел к следующему варианту.

#include <iostream>
#include <math.h>
#include <iomanip>


using namespace std;

long double Round(long double chislo, int n)
{
long double result;
__int64 iChislo;
double drob;
int k=0;
if (chislo<0)
k=-1;
else
k=1;
iChislo = (__int64) chislo;
drob = chislo - iChislo;
result=((__int64) (drob*pow(10.0,n)+k*0.5))/pow(10.0,n);
result = iChislo + result;
return result;
}


long double Cut(long double value,int m)
{
int n=0;
long double temp=abs(value);
while (temp>=1)
{
temp=temp/10;
n++;
}
value=Round(value,abs(m-n));
return value;
}

long double Sum(long double num,int n,int m)
{
long double temp=0;
long double tmp=0;
for (int i=1;i<=n;i++)
{
tmp=pow(num,i+1);
tmp=Cut(tmp,m);
for (int j=2;j<=i;j++)
{
tmp=tmp/j;
tmp=Cut(tmp,m);
}
temp+=tmp;
temp=Cut(temp,m);
}
return temp;
}

long double DeltaPogr (long double x,int n,int m)
{
long double temp=0.0;
long double an=0.0;
long double tmp=0.0;
int p=0;
int d=0;
an=pow(x,n+1);
for (int j=2;j<=n;j++)
{
an=an/j;
}
while (abs(an)>=1)
{
an/=10;
p++;
}
an=Cut(an,m);
an=an*pow(10.0,p);
cout <<" a" << n << " elem: " << an;
temp=abs(an/Sum(x,n,m));
temp=Cut(temp,m);
cout <<" Pogr iz " << n << " elem:" << temp;
tmp=abs(Sum(x,n,m));
while (tmp>=1)
{
tmp/=10;
d++;
}
tmp=Cut(tmp,m);
tmp=tmp*pow(10.0,d);
cout << " Summa : " << tmp << endl;
return temp;
}

void main()
{
long double num=0.0;
long double res=0.0;
long double x1,x2;
long double eps=0.0;
long double temp=0.0;
int m=0;
int s=0;
int N;
x1=-2.3;
x2=-17.0;
for (int i=5;i<=8;i+=3)
{
cout << "Razryadnost* manticy:" << i << endl;
N=0;
m=0;
s=0;
eps=1/pow(2.0,i-1);
temp=eps;
while (temp<=1)
{
temp*=10;
m++;
}
cout << "Epsilon:"<< eps << endl;
do
{
s++;
N++;
res=DeltaPogr(x1,s,m);
}
while (abs(res)>=eps);
cout << "Kol-vo chlenov ryada:" << N << endl;
}
}



Для точки x1 вроде правильно считает,но с выводом нечто неладное.
Volvo,Вы раньше упоминали,что выводить надо с точностью до десятой,я не совсем понял при каких случаях.Например,nxt = -4.5E-02 ,что равно -0.045 при эпсилон 0.0625 при выводе правильно выводить -0.04 или все равно выводить -0.045?Просто как я не крутился,но ответ предстает в виде.
Изображение
Если я хочу вывести все таки -0.045,то тогда у меня у всех чисел начинают появляться по 1 цифре...

Следующий вопрос,можно ли как то принудительно заставить выводить через экспоненту,потому что для второго корня числа выводятся не очень красиво.
Изображение

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

Надеюсь мне все таки ответят. smile.gif
Это текстовая версия — только основной контент. Для просмотра полной версии этой страницы, пожалуйста, нажмите сюда.