Matlab: потеря точности в расчетах. Возможно ли масштабирование переменных?

Question

Сегодня я столкнулся с проблемой точности в Matlab:

Tp = a./(3600*sqrt(g)*sqrt(K).*u.*Sd*sqrt(bB))

где

а =

                  346751.503002533 

г =

                  9.81

bB =

                  2000

Sd =

      749.158805838953
      848.621203222693
       282.57250570754
      1.69002068665559
      529.068503515487

ты =

     0.308500000000039
     0.291030000000031
      0.38996000000005
      0.99272999999926
     0.271120000000031

К =

 3.80976148470781e-009
 3.33620420353532e-009
 1.67593037457502e-008
 7.22952172629158e-005
 9.89028880679124e-009

Видимо, из-за вычислений переменных разных размеров у меня возникает проблема с точностью компьютера:

Tp =

      48.2045906674902
      48.2045906674902
      48.2045906674902
      48.2045906674902
      48.2045906674902

К сожалению, я не знаю, как справиться с этим. Я поиграл с форматом вывода, но это не проблема. Поэтому я предполагаю, что именно внутренняя точность вычислений меня привлекает. Однако, если я сам вычислю sqrt(K).* U или u.*Sd, я получу разумные значения. Только после того, как я умножу все 3 матрицы вместе, я получу одно и то же значение, хотя оно должно меняться. Я нашел этот поток, но мой случай немного отличается, потому что я не получаю произвольные значения, но они все по какой-то причине одинаковы: числовая проблема при вычислении дополнительной проекции

Я также подумал, что масштабирование всех переменных примерно так: Sd = Sd/max(Sd) может помочь, но, поскольку мне нужен очень точный и корректный по размерам результат, это не поможет.

Даже при использовании

vpa(a./(3600*sqrt(g)*sqrt(K).*u.*Sd*sqrt(bB)))

Я получаю одно и то же значение каждый раз, но с большим количеством цифр. Почему это?

Я надеюсь, что вы можете мне помочь. ура

Изменить: Вот больше кода для лучшего понимания моей проблемы:

Al = 2835000000; % [m^2]
Qp = 3000000; [m^3*s^-1]
% draw 100 uniformally distributed values for s & r
s = 600 + (8000-600).*rand(100,1);
r = 600 + (15000-600).*rand(100,1);

% calculate Sd & Rd
Sd = 680./s;
Rd = 680./r;
figure
subplot(2,1,1)
hist(Sd)
subplot(2,1,2)
hist(Rd)

%% calculate my numerically
% calculate sigma
sig = Sd./Rd;

% define starting parameters for numerical solution
t = -1*ones(size(sig));
u = zeros(size(sig));
f = zeros(size(sig));

% define step
st = 0.00001;

% define break criterion
br = -0.001;

% increase u incrementally by st until t <= br
for i=1:length(sig)
    while t(i)<br
        while t(i)<-0.1
            f(i) = sig(i)*u(i)/sqrt(pi); % calculate f for convenience
            ierfc = exp(-f(i)*f(i))/sqrt(pi) - f(i)*erfc(f(i)); % calculate integral  of complementary error function
            t(i) = (u(i)/sqrt(pi))*erfc(-f(i))*(1+sig(i))-ierfc
            u(i) = u(i) + st*1000
        end
        while t(i)>=-0.1&& t(i)<br
            f(i) = sig(i)*u(i)/sqrt(pi); % calculate f for convenience
            ierfc = exp(-f(i)*f(i))/sqrt(pi) - f(i)*erfc(f(i)); % calculate integral of complementary error function
            t(i) = (u(i)/sqrt(pi))*erfc(-f(i))*(1+sig(i))-ierfc
            u(i) = u(i) + st;
        end
    end
end
figure
hist(u)


%% calculate K from Qp
K = 3/2*pi*(Qp^(2/3)*bB^(1/3))./(g^(1/3)*u.^2.*Sd.^2*Al);

%% calculate Tp
% in hours!
Tp = (3/2*sqrt(6*pi)*sqrt(Al))./(3600*sqrt(g)*sqrt(K).*u.*Sd*sqrt(bB));

Answer 1

Я запустил этот тест, который использует только ваши векторные величины:

Sd = [                    ...
        749.158805838953  ...
        848.621203222693  ...
        282.57250570754   ...
        1.69002068665559  ...
        529.068503515487
];

u = [                     ...
        0.308500000000039 ...
        0.291030000000031 ...
        0.38996000000005  ...
        0.99272999999926  ...
        0.271120000000031 ...
];

K = [                         ...
        3.80976148470781e-009 ...
        3.33620420353532e-009 ...
        1.67593037457502e-008 ...
        7.22952172629158e-005 ...
        9.89028880679124e-009 ...
];

r = sqrt(K).*u.*Sd;
min_r = min(r);
max_r = max(r);
disp(min_r);
disp(max_r - min_r);

И я получил такой результат:

0.0143

3.2960e-17

Мне кажется, что нет реальной потери точности, но ваши векторы настроены таким образом, что они будут возвращать примерно одинаковые значения. Я имею в виду, когда значение порядка 10^−2, ошибки порядка 10^−17 довольно малы, почти что соответствуют точности представления двойных чисел (16 десятичных цифр). А двойные потери точности с плавающей запятой должны вызывать гораздо меньшую озабоченность по сравнению, например, с потерями точности при преобразовании в / из десятичного представления. Итак, вопросы: 1) надежны ли ваши источники данных и / или точны ли они? 2) вы уверены, что поэлементное произведение трех векторов не должно возвращать вектор с одинаковым значением?

LaterEdit

Мы показываем только векторную зависимость и игнорируем скаляры, потому что они вносят вклад во все компоненты вектора одним и тем же фактором. Мы будем использовать '~' для выражения пропорциональности между компонентами вектора. Тогда по вашей формуле:

K_i ~ u_i^-2 × Sd_i^-2

Tp_i ~ K_i^−1/2 × u_i⁻¹ × Sd_i⁻¹

Подключив первую формулу ко второй, вы получите:

Tp_i ~ (u_i^-2 × Sd_i^-2)^-1/2 × u_i^-1 × Sd_i^-1

или, после некоторой тривиальной алгебраической манипуляции:

Tp_i ~ u_i^{(-2× -1 / 2)} × Sd_i^{(-2× -1 / 2)} × u_i^-1 × Sd_i^-1

Tp_i ~ u_i × Sd_i × u_{i −}¹ × Sd_{i −}¹

Tp_я ~ 1_я

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