計算結果不一致の理由について

(☆)

こちらの式を計算したいです。変数はxでa,b,c,kについては定数です。

Symbolic Math toolboxで代入計算を行う方法と、Symblic Math toolboxで解析解を求めておいて後から代入計算する2つの方法を考えました。

この2つの計算方法は同じになると思ったのですが結果があわず、考えても分からないため質問させていただきます。

この理由について、分かればご教授していただけないでしょうか。

よろしくお願いいたします。

% 目的は(☆)を解くことです。

%% Symbolic Math toolboxを使って解く場合
syms x
a = 3.5e-3;
b = 2.3e3;
c = 46e-3;
k = 0.0444;

f1 = a* sinc(b*(x-c)/pi); 
f2 = f1*x;

ans= subs(int( int(f2,k,x)), k);  %(☆)
double(ans)  %  計算結果 5.7814e-09

%% Symbolic Math toolboxで解析解を求め、それに定数値を代入する
syms x a b c k
f1 = a* sinc(b*(x-c)/pi);  %sinc関数
f2 = f1*x;
sol = subs(int( int(f2,k,x)), k)
%解析解は↓のように求まる
% ((a*sin(b*(c - k)))/b + a*c*cos(b*(c - k)))/b^2 - ((c - k)*(a*cos(b*c - b*k) - a*b*c*sinint(b*(c - k)) + a*b*c*sinint(b*c - b*k)))/b^2;

% 解析解に代入
a = 3.5e-3;  
b = 2.3e3;
c = 46e-3;
k = 0.0444;
((a*sin(b*(c - k)))/b + a*c*cos(b*(c - k)))/b^2 - ((c - k)*(a*cos(b*c - b*k) - a*b*c*sinint(b*(c - k)) + a*b*c*sinint(b*c - b*k)))/b^2
%  計算結果 -2.5368e-11

Answer 1

[Hernia Baby]

おそらくですが、浮動小数点データの精度の問題だと思われます。

https://docs.microsoft.com/ja-jp/archive/blogs/jpvsblog/93

[持田]

ご回答ありがとうございます！ それも考えSymbolic(上)の方で精度を変えてみたらなにか分かるかと思いましたがよく分かりませんでした。(おそらく私の誤差に対する認識不足ですが)

Symbolic(上)の方で精度を引き上げても変わらないので、Symbolicの計算ではうまくできていると考えてもよいのでしょうか...？それともこの計算はMATLABでやるには厳しいですか？

解析解に代入する方法(下)で誤差回避する方法がもしあれば教えていただければ大変ありがたいです。

質問を重ねてしまってすみません。。 こちらの方でももう少し熟考してみます。

[Hernia Baby]

Toolboxのアルゴリズムを完全に理解してませんので誤りがあるかもしれません。

Symbolic(下)が正しい計算だと思ってます。

理由はSymbolic(上)の計算は桁の大きい状態での計算だからです。

-------------------------------------

subsのところで計算がずれています。

積分の前に数値を代入すると桁が大きすぎて精度誤差が出るのではないかと推測します。

一度分解して計算します。

Symbolic(上)の計算から

syms x
a0 = 3.5e-3;
b0 = 2.3e3;
c0 = 46e-3;
k0 = 0.0444;

f10 = a0* sinc(b0*(x-c0)/pi); 
f20 = f10*x;

sol10 =int(f20,k0,x);
sol20 =int(sol10);
sol30 = subs(sol20,x,k0);

ここで

sol30 =

(161*sinint(92/25))/50000000000 + (791*cos(92/25))/13225000000000 + (7*sin(92/25))/24334000000000

となっています。

-------------------------------------

Symbolic(下)の計算をそのあとに続けてみます

syms x a b c k
f1 = a* sinc(b*(x-c)/pi);  %sinc関数
f2 = f1*x;

% Symbolic(下)の回答に代入
sol = subs(int(int(f2,k,x)), x,k);
ans11 = subs(sol,[a,b,c,k],[a0,b0,c0,k0]);

% Symbolic(下)の回答を分解
sol11 = int(f2,k,x);
sol21 = int(sol11);
sol31 = subs(sol21,x,k);

% Symbolic(上)と同じになる
sol11 = subs(sol11,[a,b,c,k],[a0,b0,c0,k0]);
sol210 = subs(int(sol11),k0);

% Symbolic(下)と同じ
% sol21n = subs(sol21,[a,b,c,k,x],[a0,b0,c0,k0,k0]);
% sol211 = subs(sol21,[a,b,c,k],[a0,b0,c0,k0]);

% Symbolic(下)と同じ
sol31 = subs(sol31,[a,b,c,k],[a0,b0,c0,k0]);

% Symbolic(下)と同じ
sol310 = subs(sol211,k);
sol311 = subs(sol310,k0);

sol31 =

(777*cos(92/25))/26450000000000 + (7*sin(92/25))/24334000000000

です

-------------------------------------

Symbolic(上)と同じになる計算を見てみると、不定積分の前にa,b,c,kの数値を代入したときに計算が変わっています。

ここでSymbolic(上)のsol20を見てみると

sol20 =

(7*cos(2300*x - 529/5))/230000000000

- (7*sin(2300*x - 529/5))/24334000000000

+ (7*sinint(2300*x - 529/5)*(2300*x - 529/5))/230000000000

+ (7*x*sinint(92/25))/100000000

+ (7*x*cos(92/25))/10580000000

おそらくここで桁が大きすぎて誤差が出てると思われます。

-------------------------------------

結論ですが、solで出した解析解に数値を代入するのがいいと思われます。

[Hernia Baby]

symbolic math toolboxの精度について以下のリンクがありました。

16桁での計算になっているみたいです。

vpaを使うとよさそうです。

https://jp.mathworks.com/help/symbolic/increase-precision-of-numeric-calculations.html

[持田]

詳細なご返信誠にありがとうございます。

しかしながら、精度が原因ではないような気がします。

仰る通り、Symblic Math toolboxで積分した値に対して勝手に丸めが生じる現象はあると思いますが、それが根本的な原因というわけではないと思います。vpaを用いた精度引き上げは行っており、精度100桁で行いましたが、出力は変わりませんでした。(コードのPart1.参照)

また、客観性にはかけますが、こちらの計算を用いた解析を行ったところ代入が先に行われる計算結果の5.7814e-09の方が一見すると正しそうなので、私は5.7814e-09の方が正しいと思っています。(以降は、5.7814e-09が正しいと仮定して検討を行っているため前置きです)

おかげさまで、sol11からsol21の計算で差異が生じていることが明確になりました。(Part2.)

そこで、sol11からsol21の解析解を別のサイト(Wolfram|Alpha)で計算して、比較することにしました。(Part3.)

WolframAlphaによる出力→

https://www.wolframalpha.com/input/?i=Inegral%5BIntegral%5Basinc%28b%28x-c%29%29x%2C%7Bx%2Cl1%2Cx%7D%5D%2Cx%5D&lang=ja

これと、Symblic Math toolboxによって計算された解析解(sol21)を比較しますと、

後者で以下の項が余分であることが分かりました。

-(a0*c0/b0^2*cos(b0*c0-b0*k0))-a0*c0^2/b0*sinint(b0*(c0-k0)) (☆)

具体的な値を代入してみると-5.8067e-09となりまして、これはこれまでの5.7814e-09と-2.5368e-11の差分に該当します(Part.3)　つまり、原因は(☆)の項であることが分かりました。(WolframAlphaの計算結果が正しいと仮定した場合ですが)

次に、sol11からsol12の計算がSymblic Math toolbox側で間違ってしまう原因について調べるため,sol11を3つの項に分解して積分し、あとから足すことにしてみました。(Part.4)

しかしながら、この方法で行うと、5.7814e-09という正しい値が出力されました。

以上までのことをまとめると

• sol11からsol12での計算結果に誤りが生じている可能性がある
• sol11を3つの項に分けて合算するという方法であればなぜか正しい解析解がもとまる

となります。

私は予期しない数式を入力してしまっているのでしょうか。それともSymblic Math Toolboxのバグを踏んでいるのでしょうかね...途中で思い違いをしていたら申し訳ないです。

Hernia Baby はどうお考えでしょうか？もしここで解決しなければSupoprtに問い合わせることも視野にいれるつもりです。

%% Part1. 精度を引き上げても出力は変わらない
digits(100)
syms x
a = vpa(3.5e-3);
b = vpa(2.3e3);
c = vpa(46e-3);
k = vpa(0.0444);

f1 = a* sinc(b*(x-c)/pi); 
f2 = f1*x;
ans= vpa(subs(vpa(int( vpa(int(f2,k,x)))), k));  %(☆)
double(ans)  %  5.7814e-09が出力される

%% Part.2 整理
syms x a b c k
f1 = a* sinc(b*(x-c)/pi);  %sinc関数
f2 = f1*x;

% Symbolic(下)の回答を分解
sol11 = int(f2,k,x);  
sol21 = int(sol11); %ここの計算であやまりが発生している(?)

% sol11にa,b,c,kを代入してから積分
sol11 = subs(sol11,[a,b,c,k],[a0,b0,c0,k0]); 
sol210 = subs(int(sol11),k0);  
double(sol210) %  5.7814e-09

% sol21に代入
sol21n = subs(sol21,[a,b,c,k,x],[a0,b0,c0,k0,k0]); 
double(sol21n) %   -2.5368e-11

%% Part 3. symblic math toolboxでのPart.3 sol11からsol21の計算とwolframalphaでの計算との差分(原因?)
-(a0*c0/b0^2*cos(b0*c0-b0*k0))-a0*c0^2/b0*sinint(b0*(c0-k0)) % -5.8067e-09
-5.806731012715566e-09 + 5.781363279459721e-09   %-2.5368e-11

%% Part.4 sol11の解析解を積分したら誤るので、sol11の解析解を3つの項に分割して積分
syms a b c k x
f1 = (a*(cos(b*(c - k)) - cos(b*(c - x))))/b^2;
f1_int = int(f1,x);
f2 = (a*c*sinint(b*(c - k)))/b;
f2_int = int(f2,x);
f3 = - (a*c*sinint(b*(c - x)))/b;
f3_int = int(f3,x);
f_sum = f1_int + f2_int + f3_int;
sol21n = subs(f_sum,[a,b,c,k,x],[a0,b0,c0,k0,k0]); 
double(sol21n)  % 5.7814e-09が出力

[持田]

いまさらですが、最初の質問の(☆)式が間違っていますね。。

内側の積分の中にxが入ります。コードの式が合っています