0

Could anyone please tell me how to do this? I am new to Matlab as well as Mathematica. I have my mathematica coding. But, it gives different results when I run it different time. So, I want to run it in Matlab and verify my result. Please help me anyone.Really appreciate it. This contains defining functions, parametric plot and etc. I did find thru Matlab. But, I couldn't understand how to write the same program in Matlab.

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]"]]
Floris
  • 45,857
  • 6
  • 70
  • 122
TMH
  • 134
  • 1
  • 8
  • I do not get different answers when running this code repeatedly. I suspect that you have some global variable already defined. You should not assume that just because you get a "different" answer in Matlab that the Matlab answer is the right one. Numerical precision is a tricky thing, and you seem to be working in very small numbers here. – Verbeia Mar 31 '13 at 06:59

2 Answers2

0

Looks like quite a challenge for an inexperienced Matlab user. The Matlab code syntax is different from the Mathematica syntax. You should try to work on something more simple to get to know how Matlab works before trying to write such a difficult script.

Mathworks offers tutorials on their site (I've never done them, but I think it is worth the try) http://www.mathworks.nl/academia/student_center/tutorials/launchpad.html

Also, I found this comparison between the mathematica and matlab syntax, maybe it could help http://amath.colorado.edu/computing/mmm/syntaces.html

Luxus
  • 11
  • 4
  • Thank you fot the link. Meantime, I am studying Matlab. Could you please tell me how to define a function in Matlab. This is how it is in Mathematica. `SJJ[m_, a_] := 1/2 R^2 (BesselJ[m, a R]^2 - BesselJ[-1 + m, a R] BesselJ[m + 1, a R])` – TMH Mar 30 '13 at 17:07
  • function [y1,...,yN] = myfun(x1,...,xM) declares a function named myfun that accepts inputs x1,...,xM and returns outputs y1,...,yN. This declaration statement must be the first executable line of the function. Save the function code in a text file with a .m extension. The name of the file should match the name of the first function in the file. Valid function names begin with an alphabetic character, and can contain letters, numbers, or underscores. http://www.mathworks.nl/help/matlab/ref/function.html – Luxus Mar 30 '13 at 18:07
0

Using ToMatlab package you can convert Mathematica expressions to MATLAB equivalents.But it can convert things that have equivalent in MATLAB.

There is a similar question asked in Mathematica.stackexchange.

I don't think conversion of this code to MATLAB is too difficult for a beginner. MATLAB has a very useful ducumentation and you can consult it for your question asked above (about defining functions, Bessel functions, etc in MATLAB)

Community
  • 1
  • 1
Mostafa
  • 219
  • 6
  • 19
  • Thank you. I did look at the link you directed. But, when I use ToMatlab it gives me `$Failed`.Could you help me to do this? – TMH Mar 30 '13 at 18:54