0
任何人都可以告訴我如何做到這一點？我是Matlab和Mathematica的新手。我有我的mathematica編碼。但是，當我在不同的時間運行它時會產生不同的結果。所以，我想在Matlab中運行它並驗證我的結果。請幫助我的任何人。真的很感激它。 這包含定義函數，參數圖等。我確實通過Matlab找到了。但是，我不明白如何在Matlab中編寫相同的程序。在Matlab中編寫Mathematica代碼
L1 = 40;
L2 = 20;
A2 = 4.1;
D1 = 1.3;
B1 = 10;
D2 = 19.6;
B2 = 56.6;
N1 = D2 + B2;
N2 = D2 - B2;
A21 = 4.1;
F1 = (D1 \[Pi]^2)/L^2 + (B1 \[Pi]^2)/L^2;
F11 = F1 /. L -> L2;
F2 = (D1 \[Pi]^2)/L^2 - (B1 \[Pi]^2)/L^2;
\[Alpha] = D2^2 - B2^2;
\[Beta] = ((2 \[Pi]^2)/L^2 (D1 D2 - B1 B2) - 2 E1 D2 - A2^2)/\[Alpha];
\[Gamma] = (E1 (E1 - 2 (D1 \[Pi]^2)/L^2) + F1 F2)/\[Alpha];
\[CurlyPhi] = \[Pi]/180*(0);(* input angle in deg*)
\[Kappa]p2 = (-\[Beta] + Sqrt[\[Beta]^2 - 4 \[Gamma]])/2;
\[Kappa]m2 = (-\[Beta] - Sqrt[\[Beta]^2 - 4 \[Gamma]])/2;
\[Kappa]0 = Sqrt[\[Kappa]p2 /. L -> L1];
\[Kappa]01 = Sqrt[-\[Kappa]p2 /. L -> L1];
q = Sqrt[-\[Kappa]m2 /. L -> L1];
q1 = Sqrt[-\[Kappa]p2 /. L -> L2];
q2 = Sqrt[-\[Kappa]m2 /. L -> L2];
(*-----------Electron density \[DoubleStruckCapitalR](\[Rho]) \
:--------------------- *)
R\[Rho] = \[DoubleStruckCapitalN]^2 (p1*Q1*
BesselJ[m, \[Kappa]0*\[Rho]] +
p2*l1*BesselI[m, q*\[Rho]])^2 + \[DoubleStruckCapitalN]^2 (p1*
Q2 BesselJ[m + 1, \[Kappa]0*\[Rho]] +
p2* l2*BesselI[m + 1, q*\[Rho]])^2;
(*\[Rho]<R*)
R\[Rho]1 = \[DoubleStruckCapitalN]^2 (p3*\[CapitalLambda]1*
BesselK[m, q1*\[Rho]] +
p4*\[Beta]1*
BesselK[m,
q2*\[Rho]])^2 + \[DoubleStruckCapitalN]^2 \
(p3*\[CapitalLambda]2*BesselK[m + 1, q1*\[Rho]] +
p4*\[Beta]2*BesselK[m + 1, q2*\[Rho]])^2;(*\[Rho]>R*)
(*---------------Finding p1,p2,p3,p4 -----------------*)
(*
a1 p1+a2 p2+a3 p3==d1;
b1 p1 +b2 p2+b3 p3==d2;
c1 p1+c2 p2+c3 p3==d3;
Subscript[p, 1]=((Subscript[d, 3] Subscript[a, 3]-Subscript[c, 3] \
Subscript[d, 1])(Subscript[b, 2] Subscript[a, 3]-Subscript[b, 3] \
Subscript[a, 2])-(Subscript[d, 2] Subscript[a, 3]-Subscript[b, 3] \
Subscript[d, 1])(Subscript[c, 2] Subscript[a, 3]-Subscript[a, 2] \
Subscript[c, 3]))/((Subscript[c, 1] Subscript[a, 3]-Subscript[c, 3] \
Subscript[a, 1])(Subscript[b, 2] Subscript[a, 3]-Subscript[b, 3] \
Subscript[a, 2])-(Subscript[b, 1] Subscript[a, 3]-Subscript[b, 3] \
Subscript[a, 1])(Subscript[c, 2] Subscript[a, 3]-Subscript[a, 2] \
Subscript[c, 3]));
Subscript[p, 2]=(Subscript[d, 2] Subscript[a, 3]-Subscript[b, 3] \
Subscript[d, 1])/(Subscript[b, 2] Subscript[a, 3]-Subscript[b, 3] \
Subscript[a, 2])-Subscript[p, 1]((Subscript[b, 1] Subscript[a, \
3]-Subscript[b, 3] Subscript[a, 1])/(Subscript[b, 2] Subscript[a, \
3]-Subscript[b, 3] Subscript[a, 2]));
Subscript[p, 3]=Subscript[d, 1]/Subscript[a, 3]-Subscript[a, \
1]/Subscript[a, 3] Subscript[p, 1]-Subscript[a, 2]/Subscript[a, 3] \
Subscript[p, 2
];
Subscript[p, 4]=1;
*)
p1 = -((-b3 c2 d1 + b2 c3 d1 + a3 c2 d2 - a2 c3 d2 - a3 b2 d3 +
a2 b3 d3)/(
a3 b2 c1 - a2 b3 c1 - a3 b1 c2 + a1 b3 c2 + a2 b1 c3 - a1 b2 c3));
p2 = -((b3 c1 d1 - b1 c3 d1 - a3 c1 d2 + a1 c3 d2 + a3 b1 d3 -
a1 b3 d3)/(
a3 b2 c1 - a2 b3 c1 - a3 b1 c2 + a1 b3 c2 + a2 b1 c3 - a1 b2 c3));
p3 = -((-b2 c1 d1 + b1 c2 d1 + a2 c1 d2 - a1 c2 d2 - a2 b1 d3 +
a1 b2 d3)/(
a3 b2 c1 - a2 b3 c1 - a3 b1 c2 + a1 b3 c2 + a2 b1 c3 - a1 b2 c3));
p4 = 1;
a1 = \[Kappa]0 BesselJ[m, \[Kappa]0 R];
a2 = q BesselI[m, q R];
a3 = q1 BesselK[m, q1 R];
b1 = ((F1 /. L -> L1) - E1 + N1 \[Kappa]0^2) BesselJ[
m + 1, \[Kappa]0 R];
b2 = ((F1 /. L -> L1) - E1 - N1 q^2) BesselI[m + 1, q R];
b3 = -(F11 - E1 - N1 q1^2) BesselK[m + 1, q1 R];
c1 = \[Kappa]0^2 BesselJ[m + 1, \[Kappa]0 R];
c2 = -q^2 BesselI[m + 1, q R];
c3 = q1^2 BesselK[m + 1, q1 R];
d1 = -q2 BesselK[m, q2 R];
d2 = (F11 - E1 - N1 q2^2) BesselK[m + 1, q2 R];
d3 = -q2^2 BesselK[m + 1, q2 R];
(*------------Normalization constant---------------*)
\[DoubleStruckCapitalN] =
Sqrt[1/(2 Pi \
\[DoubleStruckCapitalN]1)];(*1/(\[DoubleStruckCapitalN]^2 2 Pi)=\
\[DoubleStruckCapitalN]1*)
\[DoubleStruckCapitalN]1 = (p1^2*Q1^2*SJJ[m, \[Kappa]0]) + (2*p1*Q1*
p2*l1*SJI[m, \[Kappa]0, q]) + (p2^2*l1^2*SII[m, q]) + (p1^2*Q2^2*
SJJ[m + 1, \[Kappa]0]) + (2*p1*p2*Q2*l2*
SJI[m + 1, \[Kappa]0, q]) + (p2^2*l2^2*
SII[m + 1, q]) + (p3^2*\[CapitalLambda]1^2*SKK[m, q1]) + (2*p3*
p4*\[CapitalLambda]1*\[Beta]1*
SKKab[m, q1, q2]) + (p4^2*\[Beta]1^2*
SKK[m, q2]) + (p3^2*\[CapitalLambda]2^2*SKK[m + 1, q1]) + (2*p3*
p4*\[CapitalLambda]2*\[Beta]2*
SKKab[m + 1, q1, q2]) + (p4^2*\[Beta]2^2*SKK[m + 1, q2]);
Q1 = -A2 \[Kappa]0;
Q2 = (F1 /. L -> L1) - E1 + N1 \[Kappa]0^2;
l1 = -A2 q;
l2 = (F1 /. L -> L1) - E1 - N1 q^2;
\[CapitalLambda]1 = A2 q1;
\[CapitalLambda]2 = F11 - E1 - N1 q1^2;
\[Beta]1 = A2 q2;
\[Beta]2 = F11 - E1 - N1 q2^2;
(*----------Defining the notations----------------*)
SJJ[m_, a_] :=
1/2 R^2 (BesselJ[m, a R]^2 -
BesselJ[-1 + m, a R] BesselJ[m + 1, a R])
SJJab[m_, a_, b_] :=
1/(b^2 - a^2) (a R BesselJ[m, b R ] BesselJ[m - 1, a R] -
b R BesselJ[m - 1, b R] BesselJ[m, a R])
SJI[m_, a_, b_] :=
1/(b^2 + a^2) (-a R BesselI[m, b R ] BesselJ[m - 1, a R] +
b R BesselI[m - 1, b R] BesselJ[m, a R])
SII[m_, a_] :=
1/2 R^2 (BesselI[m, a R]^2 - BesselI[m - 1, a R] BesselI[m + 1, a R])
SIIab[m_, a_, b_] :=
1/(b^2 - a^2) (-a R BesselI[m, b R ] BesselI[m - 1, a R] +
b R BesselI[m - 1, b R] BesselI[m, a R])
SKK[m_, a_] := -(1/2)
R^2 (BesselK[m, a R]^2 - BesselK[m - 1, a R] BesselK[m + 1, a R])
SKKab[m_, a_,
b_] := -(1/(
a^2 - b^2)) (b R BesselK[m, a R ] BesselK[m - 1, b R] -
a R BesselK[m - 1, a R] BesselK[m, b R])
m = 0;
R = 200;
E1 = {0.0888446, 0.153953, 0.24331};
Pm0R200 =
Plot[Piecewise[{{R\[Rho]*10^5, \[Rho] < R}, {R\[Rho]1*10^5, \[Rho] >
R}}], {\[Rho], 0, 250},
AxesLabel -> {Style["\[Rho]", Bold, FontSize -> 18],
Style["|\[CapitalPsi](\[Rho])|\!\(\*SuperscriptBox[\(\\\ \), \
\(2\)]\) (*\!\(\*SuperscriptBox[\(10\), \(-5\)]\))", Bold,
FontSize -> 15]},
BaseStyle -> {FontSize -> 15, FontWeight -> Plain,
FontFamily -> "Times New Roman"}, PlotRange -> Full,
ImageSize -> 700, PlotStyle -> Automatic,
PlotLabel ->
Style["\[DoubleStruckCapitalR](\[Rho]) vs \[Rho] : m=0 & R=200\
\[Angstrom]"]]
當我反覆運行這段代碼時，我得不到不同的答案。我懷疑你已經定義了一些全局變量。你不應該假設這只是因爲你在Matlab中得到了一個「不同的」答案，而Matlab的答案是正確的。數值精度是一件棘手的事情，你似乎在這裏工作的人很少。 – Verbeia
@Verbeia感謝您的評論。 – TMH
@Verbeia謝謝你的時間。非常感謝。 – TMH