0

今日、Matlab の精度に関する問題に遭遇しました。

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

どこ

a =

                  346751.503002533

g =

                  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

K =

 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));
4

1 に答える 1

3

ベクトル量のみを使用するこのテストを実行しました。

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 の桁のエラーはかなり小さく、double の表現精度 (10 進数で 16 桁) に近くなります。また、倍精度浮動小数点の精度損失は、たとえば、10 進数表現との間で変換する際の精度損失と比較して、はるかに懸念事項が少ないはずです。質問は次のとおりです。1) データ ソースは信頼性が高く、正確ですか。2) 3 つのベクトルの要素ごとの積が均一な値のベクトルを返すべきではないことを確信していますか?

後で編集

ベクトルの依存関係のみを示し、スカラーはすべてのベクトル コンポーネントに同じ係数で寄与するため、無視します。「~」を使用して、ベクトル コンポーネント間の比例関係を表します。次に、式に従って：

K i ~ u i −2 × Sd i −2

Tp i ~ K i −1/2 × u i −1 × Sd i −1

最初の式を 2 番目の式に代入すると、次のようになります。

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 −1 × Sd i −1

Tp i ~ 1 i

したがって、結果のベクトルTpには、すべてのコンポーネントが同じ値を持つはずです。これは、事故や精度の制限によるものではありません。これは、K、Tp、またはその両方を計算する方法によるものです。

于 2013-04-28T11:30:18.623 に答える