(物理)運動量保存はエネルギー保存則？ その17

弾性衝突に関係するエネルギーを区別し以下に列挙していきますと、
静止質量のエネルギー(mC^2)、運動質量のエネルギーがある。
運動質量のエネルギーはγmc^2のこと(特殊相対性理論では質量はγmとなることからも導出できるので。)
 運動質量のエネルギーは観察する立場により異なる相対性(並進対称性？)のエネルギーで偽？エネルギー。
弾性衝突ではこれらのエネルギーは衝突前後(というよりその時点つまり瞬間で)で保存する。エネルギー以外でも運動量、運動質量、静止質量も保存する。特殊相対性理論での近似域でも同じ。
保存量はこの三つしかない、なぜなら特殊相対性理論のエネルギーを微積分していくとエネルギー γmc^2->運動量_γmv->質量γm->運動量->エネルギーと巡回するだけだから。相対速度のBで微分してあとB->0として良い訳はこれが速度であることから例えば0の時点でも加速度は持っているからその変化は有る。
特殊相対性理論ではエネルギー、質量はc^2が有る無しだけなので等価らしい。
近似域ではエネルギーは2次式(スカラ)なのでプロットすると円、運動量は1次式(ベクトル)でプロットすると線、この交点で衝突した場合と、しなかった場合の衝突後の速度が決まる。
従って衝突前後で保存するのは当然。 もし保存しなかったらこの世界はないのだから問うこと自身意味がない。

ということで特殊相対性理論でのうエネルギーと運動量のプロットするとどうなるかみていきます。

特殊相対性理論ではエネルギーと運動量の方程式で弾性衝突後の速度は解けないので数値的に求める。
数式処理ソフトのMAXIMAで方程式の0点を飛び出させるように絶対値の逆数をとり、それを3Dグラフでプロットすれば様子がわかるだろうということです。

その様子が上の図でエネルギーの線が円で、運動量が斜め線。
結構微妙ですが、MAXIMAはそれなりに描いてくれているようです。
これをみると弾性衝突でのエネルギーと運動量の役目がよくわかります。

交わるところが解のところで、右下は衝突しなかった場合、左側が衝突した場合です。 従ってこの左を目安に数値解を求めて行きます。
MAXIMAのグラフはグリグリ動かすことができるので真上からみるとよりはっきりします。
詳しくは下に式を載せてあるのでみてください。
それで、あれやこれやを試しながら、解のあるところをある程度解の範囲に絞り込んで、それを簡易ソルバーでさらに締め上げていきます。
この簡易ソルバーは2方向から締め上げる方式で、精度は10^-4 程度が限界のようです。
四方八方から締め上げていけば精度をあげることができると思います。

下の式は順次実行しないといけないので、ctl + R で再実行するようにしてください。
maximaがおもちゃのようにヒョヒョコ再計算して行くのはみていて面白いですよ。

注意：3Dグラフをプロットするところは頭にxを付けて xplot3D( …としてMAXIMAが反応しないようにしてあるので表示したい場合はxをとってください。
----------------------------------------------------------------------------------
/* [wxMaxima: input   start ] */
kill(all);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: comment start ]
% set initial value
   [wxMaxima: comment end   ] */
/* [wxMaxima: input   start ] */
C:30;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
B:20;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
vi1:10;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
vi2:-20;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
m1:1;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
m2:2;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: comment start ]
% cov
   [wxMaxima: comment end   ] */
/* [wxMaxima: input   start ] */
vo1:(B+vi1)/(1+B*vi1/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
vo2:(B+vi2)/(1+B*vi2/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
uo1:(B+ui1)/(1+B*ui1/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
uo2:(B+ui2)/(1+B*ui2/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: comment start ]
% gamma
   [wxMaxima: comment end   ] */
/* [wxMaxima: input   start ] */
Gvi1:1/sqrt(1-vi1^2/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
Gvi2:1/sqrt(1-vi2^2/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
Gui1:1/sqrt(1-ui1^2/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
Gui2:1/sqrt(1-ui2^2/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
Gvo1:1/sqrt(1-vo1^2/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
Gvo2:1/sqrt(1-vo2^2/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
Guo1:1/sqrt(1-uo1^2/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
Guo2:1/sqrt(1-uo2^2/C^2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: comment start ]
% energy,momentum
   [wxMaxima: comment end   ] */
/* [wxMaxima: input   start ] */
Ei:Gvi1*m1*C^2+Gvi2*m2*C^2-Gui1*m1*C^2-Gui2*m2*C^2;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
Qi:Gvi1*m1*vi1*C+Gvi2*m2*vi2*C-Gui1*m1*ui1*C-Gui2*m2*ui2*C;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
Eo:Gvo1*m1*C^2+Gvo2*m2*C^2-Guo1*m1*C^2-Guo2*m2*C^2;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
Qo:Gvo1*m1*vo1*C+Gvo2*m2*vo2*C-Guo1*m1*uo1*C-Guo2*m2*uo2*C;
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(plotEi(ui1,ui2), abs(Ei));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(plotQi(ui1,ui2),abs(Qi));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(plotEo(ui1,ui2), abs(Eo));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(plotQo(ui1,ui2), abs(Qo));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(plotEQi(ui1,ui2), (abs(Ei)*abs(Qi)));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(plotEQo(ui1,ui2), (abs(Eo)*abs(Qo)));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(plotEio(ui1,ui2), (abs(Ei)*abs(Eo)));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(plotQio(ui1,ui2), (abs(Qi)*abs(Qo)));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(plotEQ(ui1,ui2),  (abs(Ei)*abs(Eo)*abs(Qi)*abs(Qo)));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: comment start ]
% plot
   [wxMaxima: comment end   ] */
/* [wxMaxima: input   start ] */
xplot3d([1/plotEQi(ui1,ui2), 0,
[ui1,-29,29], [ui2,-29,29]],[plot_format,gnuplot],[grid,50,50])$
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
xplot3d([1/plotEo(ui1,ui2), 0,
[ui1,-29,29], [ui2,-29,29]],[plot_format,gnuplot],[grid,52,52])$
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
xplot3d([1/plotEQ(ui1,ui2), 0,
[ui1,-29,29], [ui2,-29,29]],[plot_format,gnuplot],[grid,52,52])$
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
xplot3d([1/plotEQ(ui1,ui2), 500,
[ui1,-24.8,-24.4], [ui2,-0.5,0.5]],[plot_format,gnuplot],[grid,52,52])$
/* [wxMaxima: input   end   ] */
/* [wxMaxima: comment start ]
% slover
   [wxMaxima: comment end   ] */
/* [wxMaxima: input   start ] */
mysolver(slicenum,clp,cheklevel,xa,xb,ya,yb):=block(
[jj:0,kk:0,nextcheklevel,nextxa,nextya,nextxb,nextyb,dx,dy,value1,value2],
dx:(xb-xa)/slicenum, dy:(yb-ya)/slicenum,
for k:0 thru slicenum-1 do (kk:k,
    for j:0 thru slicenum-1 do (jj:j,
        nextxa:xa+jj*dx,nextya:ya+kk*dy,
        value1:float(1/plotEQ(nextxa,nextya)),
        if value1 >= cheklevel then
        (print("done left  ",jj," ",kk," ",value1," "),
        j:slicenum-1,k:slicenum-1)
  )),
for k:0 thru slicenum-1 do (kk:k,
    for j:0 thru slicenum-1 do (jj:j,
        nextxb:xb-jj*dx,nextyb:yb-kk*dy,
        value2:float(1/plotEQ(nextxb,nextyb)),
        if value2 >= cheklevel then
        (print("done right ",jj," ",kk," ",value2," "),
        j:slicenum-1,k:slicenum-1)
        )),
if value1 <= value2 then (value1:value2), nextcheklevel:value1*clp,
return([clp,value1,nextcheklevel,float(nextxa),float(nextxb),float(nextya),float(nextyb)])
)$
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
result:mysolver(100,40,1000,-24.6,-24.55,-0.2,0);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
result:mysolver(100,%[1],%[3],%[4],%[5],%[6],%[7]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
result:mysolver(100,%[1],%[3],%[4],%[5],%[6],%[7]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
result:mysolver(100,%[1],%[3],%[4],%[5],%[6],%[7]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
result:mysolver(100,%[1],%[3],%[4],%[5],%[6],%[7]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
result:mysolver(100,%[1],%[3],%[4],%[5],%[6],%[7]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
result:mysolver(100,%[1],%[3],%[4],%[5],%[6],%[7]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
result:mysolver(100,%[1],%[3],%[4],%[5],%[6],%[7]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
result:mysolver(100,%[1],%[3],%[4],%[5],%[6],%[7]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: comment start ]
check Left
   [wxMaxima: comment end   ] */
/* [wxMaxima: input   start ] */
float(plotEi(result[4],result[6]));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
float(plotQi(result[4],result[6]));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
float(plotEo(result[4],result[6]));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
float(plotQo(result[4],result[6]));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: comment start ]
check right
   [wxMaxima: comment end   ] */
/* [wxMaxima: input   start ] */
float(plotEi(result[5],result[7]));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
float(plotQi(result[5],result[7]));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
float(plotEo(result[5],result[7]));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
float(plotQo(result[5],result[7]));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: comment start ]
   [wxMaxima: comment end   ] */
/* [wxMaxima: input   start ] */
float(Gvi1*m1+Gvi2*m2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
float(Gui1*m1+Gui2*m2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(mi(ui1,ui2), float(Gui1*m1+Gui2*m2));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
mi(result[4],result[6]);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: comment start ]
   [wxMaxima: comment end   ] */
/* [wxMaxima: input   start ] */
float(Gvo1*m1+Gvo2*m2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
float(Guo1*m1+Guo2*m2);
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
define(mo(ui1,ui2), float(Guo1*m1+Guo2*m2));
/* [wxMaxima: input   end   ] */
/* [wxMaxima: input   start ] */
mo(result[4],result[6]);
/* [wxMaxima: input   end   ] */
----------------------------------------------------------------------------------