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

17 ビュー (過去 30 日間)
持田
持田 2021 年 5 月 7 日
コメント済み: 持田 2021 年 5 月 9 日
(☆)
こちらの式を計算したいです。変数は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

サインインしてこの質問に回答する。

採用された回答

Hernia Baby
Hernia Baby 2021 年 5 月 7 日
おそらくですが、浮動小数点データの精度の問題だと思われます。
https://docs.microsoft.com/ja-jp/archive/blogs/jpvsblog/93
  5 件のコメント
3 件の古いコメントを表示
持田
持田 2021 年 5 月 9 日
詳細なご返信誠にありがとうございます。
しかしながら、精度が原因ではないような気がします。
仰る通り、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が出力
持田
持田 2021 年 5 月 9 日
いまさらですが、最初の質問の(☆)式が間違っていますね。。
内側の積分の中にxが入ります。コードの式が合っています

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File Exchange で 線形代数 についてさらに検索

タグ

製品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!