умножение матриц в C и MATLAB, другой результат - PullRequest
0 голосов
/ 20 сентября 2019

Я использую 4Rungekutta для решения DGL (A x + B u = x_dot) в MATLAB и C, A - 5x5, x - 5x1, B 5x1, u 1x1, u - выводфункция синуса (2500 точек), выходные данные 4Rungekutta в MATLAB и C одинаковы до 45-й итерации, но на 45-й (в 2500 итерациях) итерации 4Rungekutta выходные данные A * x на 2-м шаге 4Rungekutta различны,матрица.я напечатал их с 30 десятичными знаками A и x одинаковы в MATLAB и C

A = [0, 0.100000000000000005551115123126,0,0,0;
-1705.367199390822406712686643004417 -13.764624913971095665488064696547 245874.405372532171895727515220642090 0.000000000000000000000000000000 902078.458362009725533425807952880859; 
0, 0, 0, 0.100000000000000005551115123126, 0;
2.811622989796986438193471258273, 0, -572.221510883482778808684088289738, -0.048911651728553134921284595293 ,0;
0, 0, -0.100000000000000005551115123126 0, 0]

x = [0.071662614269441649028635765717 ;
45.870073568955461951190955005586;
0.000002088948888569741376840423;
0.002299524406171214990085571728;
0.000098982102875767145086331744]

, но результаты A * x не совпадают, второй элемент в MATLAB равен -663.792187417201375865261070430279 C в C равен-663.792187417201489552098792046309

MATLAB
A*x = [ 4.587007356895546728026147320634
  -663.792187417201375865261070430279
  0.000229952440617121520692600622
  0.200180438762844026268084007825
  -0.000000208894888856974158859866];

C
A*x = [4.587007356895546728026147320634
 -663.792187417201489552098792046309
  0.000229952440617121520692600622
  0.200180438762844026268084007825
  -0.000000208894888856974158859866];

хотя разница небольшая, но мне нужен этот результат, чтобы получить конечную разницу, в этот момент результат будет более очевидным

кто-нибудь знает почему?

1 Ответ

1 голос
/ 22 сентября 2019

Сколько цифр вы считаете нужным?У вас одинаковые первые 16 цифр каждого числа равны, что является приблизительным количеством данных, которые double обычно могут представлять и хранить внутри.Вы не можете получить больше, даже если вы заставите свои процедуры печати печатать больше цифр, они будут печатать мусор.То, что происходит, это то, что вы сказали, чтобы получить, скажем, 120 цифр на ваши процедуры печати ... и они будут печатать их, как правило, умножая остаток (что бы это ни было). Поскольку числа представлены в базе 2, вы обычно не получаетенули однажды прошли внутреннюю прецизионность числа ... и реализациям печати не нужно согласовывать напечатанные цифры, если у вас больше цифр не представлено в вашем номере.

Предположим, вы на мгновениеиметь ручной калькулятор, который имеет только 10 цифр точности.И вам даны цифры из 120 цифр.Вы начинаете вычислять и получаете результаты только с 10 цифрами ... но вас попросили напечатать отчет с результатами из 120 цифр.Ну .... так как общий расчет не может быть сделан с более чем 10 цифрами, что вы можете сделать?Вы используете калькулятор, который не может дать вам требуемое количество цифр ... и более того, число из 10 базовых цифр в 52-битном значении не является целым числом цифр (в 52-битном значении есть 15.65355977452702215111442252567364 десятичных цифр и),Что вы можете сделать, вы можете заполнить нулями (неверно, скорее всего), вы можете заполнить эти места рубинами (что никогда не повлияет на конечный 10-значный результат) или вы можете пойти в Radio Shack и купить 120-значный калькулятор.Процедуры печати с плавающей запятой используют счетчик, чтобы указать, сколько раз зайти в цикл и получить другую цифру, они обычно останавливаются, когда счетчик достигает своего предела, но не прилагают никаких дополнительных усилий, чтобы узнать, сошли ли вы с ума и указалибольшое количество цифр ... если вы попросите 600 цифр, вы просто получите 600 итераций цикла, но цифры будут поддельными.

Вы должны ожидать разницу в одной части в 2^52 в doubleчисло, так как это число двоичных цифр, используемых для значимого (это приблизительно 1009 *, поэтому вам нужно умножить это число на число, которое вы вывели, чтобы увидеть, насколько велика ошибка округления, приблизительно), если вы умножаетечисло 663.792187417201375865261070430279 тем самым вы получите 1.473914740073748177152126604805902e-13, который является оценкой того, где в числе находится последняя действительная цифра этого числа.Вероятно, оценка ошибки будет намного больше из-за большого количества умножений и сумм, необходимых для вычисления ячейки.В любом случае разрешение 1.0e-13 очень хорошее (субатомное различие, если значения будут длинами и единицами в метрах).

РЕДАКТИРОВАТЬ

в качестве примера, просто рассмотрите следующую программу:

#include <stdio.h>
int main()
{
        printf("%.156f\n", 0.1);
}

если вы запустите его, вы получите:

0.100000000000000005551115123125782702118158340454101562500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

, что действительно является наиболее (точным) приближением к внутреннему представлению числа 0.1, которое может представлять машинав базе 2 с плавающей запятой (0.1 является периодическим числом, если оно представлено в базе 2). Его представление:

0.0001100110011(0011)*

, поэтому оно не может быть точно представлено 52 битами, повторяющими шаблон 1100на неопределенный срок.В какой-то момент вы должны сократить, и подпрограмма printf продолжает добавлять нули вправо, пока не дойдет до представленного выше представления (все конечные цифры в базе 2 представляются в виде конечного числа цифр в базе 10, нообратное неверно (потому что все факторы 2 в 10, но не все факторы 10 в 2).

Если вы поделите разницу между 0.1 и 0.1000000000000000055511151231257827021181583404541015625 между 0.1 вы получите 5.55111512312578270211815834045414e-17, что составляет примерно 1/2 ^ 54 или одну четверть (примерно одну пятую) от лимита 1/2^52, который я показал вам выше. Это самое близкое число, представимое с 52 битами от числа 0.1

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...