数式を積分する時に発生するエラー

物理現象を数式で表そうとしているのですが、

途中の積分でエラーが発生します。

物理現象はすべて真空中で発生するものであり、

タングステン線に電気を流してそこから電子を取り出し、これに指向性を持たせて、真空中に安置された絶縁体の粉体に照射します。

この時、安置された粒子の表面には隣接する粒子同士によって構成される小さな空洞（微小空洞と呼びます）が存在しており、

そこに電子が照射されると、微小空洞内でさらに電子が発生するというものです。

以下にコードを示します。

（前略）

%T=[100 300 700 1200]+273.15%aに対応したタングステンフィラメント温度
%a=[7.3 12.4 24 39]%体積抵抗率_[μΩcm]_八光電機ホームページより
T=700
a=24

l=0.01%タングステン長さ[m]
D=0.1%タングステン直径[mm]
S=(0.1/2)^2*pi%断面積_mm2
R=a*(l/S)*0.01%抵抗値_Ω_八光電機ホームページより_この計算において、単位はこのままでよいらしい。
%I=[1;3;5;7;10]%タングステンにかかる電流
I=3%タングステンにかかる電流
V=I*R%タングステンにかかる電圧

A=120*10^6%熱電子放出係数_[A・m^-2K^-2]
e=1.602176634*10^-19%電機素量_C
Vfai=4.52%仕事関数_eV
kB=1.3806488*10^-23%ボルツマン定数_J・K/mol
J1=A.*(T.^2).*exp(-(e*Vfai)./(kB.*T))%熱電子放出電流密度_A/m2

J1enum=J1/e;%J1中の電子数_個
Vg1=40;%電子銃のグリッド電圧_V
eV=J1enum*Vg1;%電子ビームのエネルギー_eV
Te=(2*eV)/(3*kB);%電子ビームの電子温度_eV

ep=8.85418781*10^-12;%真空の誘電率
Tee=2.5;%二次電子温度;
r=1.0*10^-6;%ダスト粒子の半径_m
s=2*r;%微小空洞の直径的大きさ
S=0.88%表面清浄度0.75~0.88 [-]
Cf=(4.3*10^-20)/0.0000000000000000033;%接触定数_指定された数
E0=0;%初期電場
et=(s/Tee)*sqrt(((2*Cf*(S^2))/(r*pi*ep)-2*E0^2));%二次電子の放電電荷数η
fai=(et*Tee)/e;%粒子の表面電位

syms t
Em=450;%一次電子の最大降伏エネルギー_450eV(数値は指定されたもの。)
x=@(t) (t^5)*exp(-(Em/(4*kB*Te))*(t^2)+t);
 Jpart4=(Em./(4*kB*Te)).^2
% G1=exp(-2.034678788509273e+59)　計算の確認用
%  G2=exp-9.082219408774468e+15　計算の確認用
%  G3=exp-2.319180660587761e-18　計算の確認用
%  G4=exp-5.689743803633560e-35　計算の確認用
 
G=exp(-Em./(4*kB*Te))　%デバッグ用
H=G*(t.^2)+t　　　　　　%デバッグ用

sekibunx=integral(x,inf,0);←現在のエラー発生個所
delm=1.2;%最大二次δm_指定されたもの_(Halekas et as., 2008)

Jpart1=3.7*delm*J1
Jpart2=-abs(e*fai)/(kB*Te)
Jpart3=exp(Jpart2)
%expは、関数電卓ではe^230程度以上ではエラー、e^-220以下では0となる。
Jpart4=(Em/(4*kB*Te))^2
%Jpart_all=3.7*delm*J1*exp(-abs(e*fai)/(kB*Te))*(Em/(4*kB*Te))^2　

%以下、Jpart_all=Jsecです。計算を慎重に行いたかったためにJpar1~Jpart4,Jpart_allを設けました。
Jpart_all=Jpart1*Jpart3*Jpart4%二次電子電流密度_A/m2
%J_sec=3.7*delm*J1*exp(-abs(e*fai)/(kB*Te))*(Em/(4*kB*Te))^2*sekibunx;%二次電子電流密度_A/m2
Jsec=Jpart_all*sekibunx
%J2=%微小空洞の帯電時間を表す式にぶち込む数式

このコードを走らせた時、以下のようなエラーが出ます。

エラー:  ^  (line 51)　%以下のエラーで表示されている(line 〇〇)は、上記のコードで前略とされている分ズレます。
行列をべき乗にするには次元が正しくありません。行列が正方行列で、べき指数がスカラーであることを確認
してください。行列を要素ごとにべき乗するには、'.^' を使用してください。

エラー: effect_of_beam_current>@(t)(t^5)*exp(-(Em/(4*kB*Te))*(t^2)+t)

エラー: integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

エラー: integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

エラー: integralCalc (line 83)
        [q,errbnd] = vadapt(@AToInfInvTransform,interval);

エラー: integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

エラー: effect_of_beam_current (line 66)
sekibunx=integral(x,inf,0);

自分で試したこと

実は現在のエラーが発生する前、計算は何かがおかしく正しい計算がされない状態であったもののコード自体は全部走らせることはできました。

コードを拡張してTとIを行列にしての計算を試みました。

結果、同じ場場所で次元が違うとか正方行列になっていないのではとか出てきたので、

仕方なく現在のように行列として表記しないようにしています。

しかしながら現在も「行列をべき乗にするには次元が正しくありません。行列が正方行列で、べき指数がスカラーであることを確認

してください。行列を要素ごとにべき乗するには、'.^' を使用してください。」

と、行列を使用していないにも関わらず表示されているため非常に困惑しています。

不明な点がありましたら対応致します。

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

Answer 1

[stozaki]

Noruji Muto様

エラーメッセージの以下の箇所に着目してみました。

エラー: effect_of_beam_current>@(t)(t^5)*exp(-(Em/(4*kB*Te))*(t^2)+t)

先ず

T = 700
.
.
.
fai=(et*Tee)/e;%粒子の表面電位

syms t
Em=450;%一次電子の最大降伏エネルギー_450eV(数値は指定されたもの。)
x=@(t) (t^5)*exp(-(Em/(4*kB*Te))*(t^2)+t);

上記 fai までの定義が正しいとして、関数ハンドル "x" を作成している式の "演算子" に問題があるように見えます。

引数 t の演算に関連する演算子は行列演算ではなく、要素毎の演算となるのではないでしょうか？ (以下のように)

x=@(t) (t.^5).*exp(-(Em/(4*kB*Te))*(t.^2)+t);

上記で発生していたエラーは取り除くことができます。

stozaki

[Noruji Muto]

ありがとうございます。

確認できました！