今日、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));