Manipulate[ Plot[a*Sin[n*x], {x, 0, 2*Pi}], {{n, 1}, 1, 20, Appearance -> "Labeled"}, {{a, 1}, -10, 10, Appearance -> "Labeled"} ]
(* ちょっと複雑な例 *) inSin[A_, w_, k_, t_, x_] := A*Sin[w*t + k*x] outSin[A_, w_, k_, t_, x_] := A*Sin[w*t - k*x + Pi] cwInSin[\[Rho]0_, A_, w_, k_, t_, x_] := \[Rho]0/(1 + k*A*Cos[w*t + k*x]) cwOutSin[\[Rho]0_, A_, w_, k_, t_, x_] := \[Rho]0/(1 - k*A*Cos[w*t - k*x + Pi]) cwInSinPlusOutSin[\[Rho]0_, A_, w_, k_, t_, x_] := \[Rho]0/(1 + k*A*Cos[w*t + k*x] - k*A*Cos[w*t - k*x + Pi]) Manipulate[ Plot[{ inSin[A, w, k, t, x], outSin[A, w, k, t, x], inSin[A, w, k, t, x] + outSin[A, w, k, t, x], x, (cwInSin[\[Rho]0, A, w, k, t, x] - \[Rho]0), (cwOutSin[\[Rho]0, A, w, k, t, x] - \[Rho]0), (cwInSinPlusOutSin[\[Rho]0, A, w, k, t, x] - \[Rho]0), (\[Rho]0) }, {x, -3, 12}, PlotRange -> {-0.74*range, 0.74*range}, AspectRatio -> ratio, ImageSize -> 600, PlotStyle -> { {Red, Thickness[0.0003]}, {Green, Thickness[0.0003]}, {Blue, Thickness[0.003]}, {Black, Thickness[0.0003]}, {Magenta, Dashed, Thickness[0.002]}, {RGBColor[0, 0.7, 0], Dashed, Thickness[0.002]}, {Cyan, Thickness[0.004]}, {Black, Dashed, Thickness[0.0003]} }, PlotLegends -> { "inSin 入射波", "outSin 反射波", "inSin+outSin 合成波", "x 参考", "cw InSin 入射波の密度変化", "cw OutSin 反射波の密度変化", "cw InSin+OutSin 合成波の密度変化", "\[Rho]0 波が無いときの媒質密度" }], {{ratio, 0.5, "ranito(縦横比を変更)"}, 1/10, 1, Appearance -> "Labeled"}, {{\[Rho]0, 3, "媒質の密度"}, 0, 5, Appearance -> "Labeled"}, {{range, 5, "range(縦軸幅を変更)"}, 1, 20, Appearance -> "Labeled"}, {{A, 0.45, "振幅(入射波,反射波)"}, 0, 1/(2*k), Appearance -> "Labeled"}, {{w, 1, "角振動数"}, 0, 10, Appearance -> "Labeled"}, {{k, 1, "波数"}, 0, 2, Appearance -> "Labeled"}, {{t, 0.544, "時刻"}, 0, 2*Pi/w, Appearance -> "Labeled"}, ControlPlacement -> {Top} ]
V[x_, y_, z_] := -1/Sqrt[x^2 + y^2 + z^2] (* クーロン引力 *) m = {V[x, y, z], D[V[x, y, z], x], D[D[V[x, y, z], x], y], D[D[D[V[x, y, z], x], y], z]} (* 3階微分までの結果を、行列mに代入 *) mm = Simplify[m /. {y -> 0, z -> 0}, Assumptions -> {x > 0}] (* 正のx軸上の結果を、行列mmに代入 *) {mm[[2]], mm[[1]]} (* 確認のため∂xVとVを行列要素要素として表示 *) mm[[2]]/mm[[1]] (* ∂xV / V の比 *)
$Context (*現在のコンテキスト(名前空間のようなもの)を表示する*) Global` (*デフォルトで使われているコンテキストらしい*) 変数等の正式名称は、 コンテキスト名`簡略名 である。現在のコンテキスト内の変数等であれば簡略名だけで使用可。
Manipulate[ Play[Sin[2*Pi*(440 + s)*t], {t, 0, 1}], {{s, 0}, -100, 100, Appearance -> "Labeled"}]
DSolve[{y''[x] + \[Lambda]* y[x] == 0, y[0] == 0, y[\[Pi]] == 0}, y[x], x] (*二階微分方程式を解いた例*)
帯域的な仮定を設定する例) $Assumptions = {x \[Element] Reals && y \[Element] Reals && z \[Element] Reals && Rx \[Element] Reals && Ry \[Element] Reals && Rz \[Element] Reals && A > 0 && a > 0 && b > 0 && r >= 0}
ParametricPlot[{ {Re[ArcTan[4*Exp[I*\[Theta]]]], Im[ArcTan[4*Exp[I*\[Theta]]]]}, {Re[ArcTan[2*Exp[I*\[Theta]]]], Im[ArcTan[2*Exp[I*\[Theta]]]]}, {Re[ArcTan[1*Exp[I*\[Theta]]]], Im[ArcTan[1*Exp[I*\[Theta]]]]}, {Re[ArcTan[1/2*Exp[I*\[Theta]]]], Im[ArcTan[1/2*Exp[I*\[Theta]]]]}, {Re[ArcTan[1/4*Exp[I*\[Theta]]]], Im[ArcTan[1/4*Exp[I*\[Theta]]]]}, {Re[Exp[I*\[Theta]]], Im[Exp[I*\[Theta]]]}}, {\[Theta], 0, Pi*2}, PlotLegends -> "Expressions", PlotRange -> {{-2, 2}, {-2, 2}}] ParametricPlot3D[{Re[ArcTan[r*Exp[I*\[Theta]]]], Im[ArcTan[r*Exp[I*\[Theta]]]], Abs[ArcTan[r*Exp[I*\[Theta]]]]}, {\[Theta], 0, Pi*2}, {r, 0, 4}, AxesLabel -> {"Re", "Im", "Abs"}, PlotRange -> {{-2, 2}, {-2, 2}, {0, 2}}, Mesh -> Automatic] Manipulate[ ParametricPlot[{{Re[ArcTan[r*Exp[I*\[Theta]]]], Im[ArcTan[r*Exp[I*\[Theta]]]]}, {Re[r*Exp[I*\[Theta]]], Im[r*Exp[I*\[Theta]]]}}, {\[Theta], Pi/2, Pi*3/2}, PlotLegends -> {"ArcTan[r*Exp[I*\[Theta]]]", "r*Exp[I*\[Theta]]"}, PlotRange -> {{-2, 2}, {-2, 2}}], {{r, 0.5}, 0, 10, Appearance -> "Labeled"}]
ver.11では次のように代替できそうである。
(* 複素関数のプロット 2021/12/13 *) Options[ParametricPlot3DComplex] = { "inc" -> 0, (* Arg[]の結果にincを足す。z^(1/n)などの単純な場合のみ、リーマン面の葉を描くのに有効だろう。0は主分枝を描く。 *) "pRange" -> {Automatic, Automatic, {0, 5}}, (* プロットする範囲{実軸,虚軸,Abs軸} *) "LogAbs" -> "LogAbs", (* プロットの縦軸を対数にするかなど, Abs, LogAbs, LogLogAbs *) "domainScale" -> 1(* 定義域をdomainScale倍に広げる *) }; ParametricPlot3DComplex[funcz_, OptionsPattern[]] := Module[{myAbs, myAxesLabel}, Switch[ OptionValue["LogAbs"], "LogLogAbs",{(* Log Log scale *) myAbs = Log[Log[Abs[funcz /. {z -> x + I*y}]]]; myAxesLabel = {"Re", "Im", "LogLogAbs"};}, "LogAbs", {(* Log scale *) myAbs = Log[Abs[funcz /. {z -> x + I*y}]]; myAxesLabel = {"Re", "Im", "LogAbs"}; }, _, {(* Normal scale *) myAbs = Abs[funcz /. {z -> x + I*y}]; myAxesLabel = {"Re", "Im", "Abs"}; } ]; ParametricPlot3D[{x, y, myAbs}, {x, -3*OptionValue["domainScale"], 3*OptionValue["domainScale"]}, {y, -3*OptionValue["domainScale"], 3*OptionValue["domainScale"]}, PlotRange -> OptionValue["pRange"], AxesLabel -> myAxesLabel, ViewPoint -> {0, 0, Infinity} (*Above*), ColorFunction -> Function[{px, py, pz}, Hue[(Arg[funcz /. {z -> px + I*py}] + OptionValue["inc"])/(2*Pi)]], ColorFunctionScaling -> False, MeshFunctions -> {Function[{px, py, pz}, (Arg[funcz /. {z -> px + I*py}] + OptionValue["inc"])/(2*Pi)], Function[{px, py, pz}, pz]}, PlotStyle -> Opacity[1.0], Mesh -> 10, PlotPoints -> 30, MaxRecursion -> 5, Lighting -> {{"Ambient", RGBColor[.8, .8, .8]}}]] (* 複素関数のプロットのカラーマップ *) ParametricPlot3DComplexColormap[] := Module[{anglelist}, anglelist = {1, 1 + I, I, -1 + I, -1, -1 - I, -I, 1 - I, 1 - I/1000}; MatrixForm[{ Table[Hue[N[Arg[e]/(2*Pi)]], {e, anglelist}],(*Hue[負値]は、負値に整数が加えられ0~1に範囲になるようだ*) Table[N[Arg[e]/(2*Pi)], {e, anglelist}], Table[Arg[e], {e, anglelist}]}]]
実行例
ParametricPlot3DComplexColormap[]
ParametricPlot3DComplex[Log[z], "inc"->0, "pRange"->{Automatic, Automatic, {0, 5}}, "LogAbs"->"Abs", "domainScale"->1] や ParametricPlot3DComplex[Log[z]]
ParametricPlot3DComplex[Exp[z], "inc"->0, "pRange"->{Automatic, Automatic, {0, 7}}, "LogAbs"->"Abs", "domainScale"->2.09]
(注)黒い波線は数値計算上の問題だろう
ParametricPlot3DComplex[z, "inc"->0, "pRange"->{Automatic, Automatic, {0, 5}}, "LogAbs"->"Abs", "domainScale"->1]
[c.f.] https://mathrelish.com/mathematics/singularity-in-complex-analysis
Function[{x, y}, 10*x + y][a, b] または (10*#1 + #2)&[a, b]
#サイズの異なる演算はブロードキャストされる(pythonと同じ) {1, 2} + 3*f[{1, 2}] #{1, 2} + 3*( {1, 2}+{100, 200} )の計算 4 + 3*f[4] #{4, 4} + 3*( {4, 4}+{100, 200} )の計算
Limit[Abs[x - Rx]/(x - Rx), x -> Rx] Limit[Abs[x - Rx]/(x - Rx), x -> Rx, Direction -> "FromAbove"] #上から Limit[Abs[x - Rx]/(x - Rx), x -> Rx, Direction -> "FromBelow"] #下から
ある数式aがあるとして、Sqrt[X^2+y^2]をrに置き換えたいとき、 a /. {Sqrt[x^2 + y^2] -> r, 1/Sqrt[x^2 + y^2] -> 1/r} とする(a内に 1/Sqrt[X^2+y^2] があっても、ひとつめの規則で1/rに置き換えられないため)。
psi[x0_, px0_, w_, x_, t_] := 1/((I*t/w + 2*w)^2*Pi/2)^(1/4) * Exp[(t/(2*w^2) + I)/(t^2/w^2 + 4*w^2)*I*(x - x0 - px0*t)^2 + I*px0*(x - x0 - px0/2*t)] rho[x0_, px0_, w_, x_, t_] := 1/Sqrt[Pi/2*t^2/w^2 + 2*Pi*w^2] * Exp[-2/(t^2/w^2 + 4*w^2)*(x - x0 - px0*t)^2] Manipulate[ Plot[ { rho[2, px0, w, x, t], Re[psi[2, px0, w, x, t]], Im[psi[2, px0, w, x, t]] }, {x, 0 - 1, 15 + 1}, ImageSize -> 700, GridLines -> Automatic, PlotLegends -> "Expressions", PlotRange -> {-0.8, 0.8} ], { {px0, 2.1602464}, 0, 5, Appearance -> "Labeled"}, { {w, 0.5}, 0, 1.5, Appearance -> "Labeled"}, { {t, 0}, 0, 5, Appearance -> "Labeled"} ]
wd[k_] := k^2/2 (*wd[k_]:=k*) wav1[k_, w_, x_, t_, e_] := Im[Exp[I*(k*x - w*t + e)]] wav2[k_, w_, x_, t_, e_] := Im[Exp[I*(k*x - w*t + e)]] Manipulate[ Plot[{ wav1[k1, wd[k1], x, t, 0], wav2[k2, wd[k2], x, t, e2], wav1[k1, wd[k1], x, t, 0] + wav2[k2, wd[k2], x, t, e2] }, {x, -50, 50}, PlotRange -> {-5, 5}, ImageSize -> 600, PlotStyle -> {{Red, Thickness[0.0003]}, {Blue, Thickness[0.0003]}, {Green, Thickness[0.003]}}, PlotLegends -> {"wav1", "wav2", "wav1+wav2"}], {{k1, 1.0}, -5, 5, Appearance -> "Labeled"}, {{k2, 0.8}, -5, 5, Appearance -> "Labeled"}, {{e2, 0.0}, -Pi, Pi, Appearance -> "Labeled"}, {{t, 0}, 0, 100, Appearance -> "Labeled"} ]
「一般的な場合」 ファイル→印刷設定→印刷オプションで、 セルのブラケットを印刷と、領域マークを印刷のチェックをする、 マージン設定で上下左右の余白を5mmにする(左右の余白が大きすぎるため) →A3でpdfとして印刷 →AcrobatでA4に縮小印刷
「長い数式を印刷する場合」 ファイル→印刷設定→印刷用環境 →Condensedを選ぶ(表示画面と似た感じで、やや小さめのカラーで印刷できる) (Printoutだと白黒で行間がもう少し狭めになる) ファイル→印刷設定→印刷オプション →マージンを 20mmに小さくする。 A2サイズ横で、pdf出力してから、pdfをA4サイズに縮小して印刷する。
虚数単位 I ネイピア数 E 円周率 Pi 無限大 Infinity
Exp, Log Sqrt Abs, Arg Sign Sin, Cos, Tan
\[Alpha] で α になる \[Theta] で θ
下付き文字をつけたい文字を選択か文字の直後で、Ctrl+マイナス(-) を押す→プレースフォルダ があらわれる
上付き文字をつけたい文字(べき乗になる)を選択し、Ctrl+^ を押す →プレースフォルダ
分子の入力後、Ctrl+/ で分数の分母を二次元式で入力できる、らしい
Ctrl-Enter で列ベクトルの次の行の入力ができる
InputForm["hogehoge"] # 数式hogehogeを、一次元の形式(通常の入力に使われる形式)で表示 OutputForm["hogehoge"] # 通常の画面に出る形式で表示 StandardForm["hogehoge"] # see manual(たぶん通常はOutputFormと同じ) FullForm["hogehoge"] # see manual(たぶん使うことはない)
Style[ hogehoge, Smaller ] やや小さいフォントで出力する
N[ hogehoge, [桁数] ]//InputForm # hogehogeを数値として表示する。 # //InputFormを付けると一次元の入力形式で得られるが、 # 指定した桁数以上(25~50桁など)が出力される。 # 10のべき乗部分は、`桁数*^10 となる。
Round[...] #数値で、四捨五入して出力 Floor[...] #数値で、切り捨てして出力 Ceiling[...] #数値で、切上げして出力
Limit[expr, x->a] #exprのx->aの極限を計算する
Series[ expr, {x,x0,n} ] #x=x0でn次までexprをテイラー展開する Normal[ Series[ expr, {x,x0,n} ] ] #余剰項を除く形を得る
D[] 微分。高階微分はmanualをみよ
Integrate[ ] 積分。重積分はmanualをみよ 仮定付き積分の例) Integrate[s*Exp[-s^2/2], {s, 0, x}, Assumptions-> {x > 0}]
Solve[ 式, 解く変数 ] x /. Solve[ 式, x ] (*xの解のみのリスト出力する (Solveは変換規則のリストとして解を出力する。そのため、余分な所を消すのに置換表示している)*)
例) Part[z /. Solve[1 - 3*z^2 + 2*z^3 == x, z], 1 ;; 3] (*方程式をzについて解き、解のみを1~3個出力する*)
?hoge # hogeの説明を表示 ?*hoge # *hogeの関数名一覧を表示(ワイルドカードも利用可)
Options[hoge] # hoge関数で利用できるオプションの一覧表示
hoge; #式hogeの末尾に、セミコロンを付ける
/. 置換表示 -> 変数規則 例) x^2 /. x -> 2
assump = Assumptions -> {x0 ∈ Reals && px0 ∈ Reals && t >= 0 && w > 0 && x ∈ Reals && k ∈ Reals } などする。そして、 Integrate[ func[x], {x,-Infinity,Infinity}, assump ] Simplify[ func[x], assump ] で、仮定を利用する。
ClearAll["Global`*"] #何もかもをクリアする ClearAll[x] #変数xの何もかもをクリアする 詳しくはmanualをみよ。 Clear[ hoge ] 値、定義をクリア Clear["Global`*"] 現行のセッションで作成された値、定義をすべてクリア ClearAll[ hoge ] 記号の値、定義、属性、メッセージ、デフォルト設定をクリア Remove[ hoge ] シンボルの名前が認識されないように、完全に除去
<< Utilities`CleanSlate` #この2行で、Clear["Global`*"]をし、In/Out番号も初期化し、 CleanSlate[] #開放メモリ量のレポートも出すらしい ?Utilities`CleanSlate`* #詳細はこちら
m={ {1,2,3}, {4,5,6} }
Table[ x^i*y^j, {i, 1, 3}, {j, 1, 2}]
m=DiagonalMatrix[{1,2,3,4}] (* 対角行列 *)
m[[1,All]] (* mの1行目だけ取り出す。mは2次元配列より大くても同様に指定可 *) Part[m,1,All] (* 〃 *)
m[[1;;1, 1;;2]] (* 1-1行目x1-2列目からなる部分行列 *)
m // MatrixForm (* 行列表記で出力 *) m // TableForm (* テーブルで出力,要素を並べる *) または MatrixForm[m] TableForm[m] なお、//MatrixForm をつけた形で行列を、変数に代入すると、 その後の行列積などが計算できない。
A.B (* 積 *) Transpose[A] (* 転置 *) Inverse[A] (* 逆行列 *) MatrixPower[A,n] (* 行列Aの n乗 *) A*B (* 行列A,Bの各要素をお互いに掛ける *) A^n (* 行列Aの各要素を n乗する *)
ベクトルは、1行または1列の行列として扱う。
Norm[v] ベクトルのノルム、行列のノルム(詳細は調べよ)
(*例1*) Plot[{ b*(Exp[k*(b - a)*t] - 1)/(b/a*Exp[k*(b - a)*t] - 1) /. { a -> 1, b -> 0.5, k -> 1}, b*(Exp[k*(b - a)*t] - 1)/(b/a*Exp[k*(b - a)*t] - 1) /. { a -> 1, b -> 2, k -> 1}, a*a*k*t/(a*k*t + 1) /. { a -> 1, k -> 1} }, {t, 0, 10} ]
(*例2*) omega1 = 0.00288; E01 = 7.76; (*V/Ang*) Eenv[t_] := Which[ t < 2 Pi / omega1, omega1*t/2/Pi, t < 4 Pi / omega1, 2 - omega1*t/2/Pi, True, 0] Elaser1[t_] := Eenv[t]*(E01*Sin[omega1*t]) Plot[{Elaser1[t], Eenv[t], E01*Sin[omega1*t]}, {t, 0, 2500*100/10}, PlotLegends -> "Expressions", PlotRange -> {-10, 10}, GridLines -> Automatic, AxesLabel -> {"au(time)", "V/Ang"}, PlotStyle -> {Thickness[0.01], Thickness[0.01], {Dashed, Thickness[0.0001]}}]
(*例3*) omega2 = 0.00144; E02 = 2.74 ;(*V/Ang*) Elaser2[t_] := (Sin[omega2*t/16])^2*(E02*Cos[omega2*t]) Plot[{Elaser2[t], (Sin[omega2*t/16])^2, E02*Cos[omega2*t]}, {t, 0, 2500*100/3}, PlotLegends -> "Expressions", PlotRange -> {-3, 3}, GridLines -> Automatic, AxesLabel -> {"au(time)", "V/Ang"}, PlotStyle -> {Thickness[0.01], Thickness[0.01], {Dashed, Thickness[0.0001]}}]
PlotStyle -> { Thickness[0.005] } #0.005がデフォルトに近い
PlotLegends -> "Expressions"
PlotStyle -> { Thickness[0.01], Thickness[0.01], {Dashed, Thickness[0.005]}, {Dashed, Thickness[0.005]} } #4本分指定
PlotStyle -> { {Blue, Thickness[0.01]}, {Orange, Thickness[0.01]}, {Blue, Dashed, Thickness[0.005]}, {Orange, Dashed, Thickness[0.005]}}
例えば、 Plot[ {0.3*Sin[x],x}, {x, -3, 12}, PlotRange -> {-0.74, 0.74}, AspectRatio -> 1/10, ImageSize -> Full ] # 横軸幅 15、縦軸幅 1.5、縦横比が 1.5/15=1/10なので、スケールは同じくなる。 # ImageSize->Fullを付けるとウィンドウ幅に合わせて拡大されたグラフが描かれる。
(* z^(1/3)=Exp[1/3 Log[Abs[z]] + I/3 Arg[z]+I/3* 2 m Pi where -Pi<=Arg[z]<Pi and m=0,1,2 のリーマン面を描く *) Plot3D[ { Arg[( x + I*y)^(1/3)] + 0* 2 Pi /3, (*m=0,mathematicaでの主値の範囲*) Arg[( x + I*y)^(1/3)] + 1* 2 Pi /3, (*m=1*) Arg[( x + I*y)^(1/3)] + 2* 2 Pi /3 (*m=2*) }, {x, -2, 2}, {y, -2, 2}, AxesLabel -> Automatic, (*軸ラベルを表示*) PlotLegends -> "Expressions" (*凡例表示。ひとつだと出ない*), LabelStyle -> Directive[Bold, Medium], (*軸ラベルを太く大きく*) Mesh -> 3, PlotStyle -> {Opacity[0.8]}, (*透過度。色も可*) ExclusionsStyle -> {None, Red}, (*値が飛ぶ所を描かない・端を赤色にする*) Ticks -> {Automatic, Automatic, {-Pi, -Pi/2, 0, Pi/2, Pi, 3/2 Pi}}, (*軸目盛*) ViewPoint -> {-5, -3.4, 3}*1000 ]
z[r_, t_] := r Exp[ I t ]
S[n_, x_] := Which[ n == 0, Which[ x < -1, 0, x < 0, 1 - 3*x^2 - 2*x^3, x < 1, 1 - 3*x^2 + 2*x^3, True, 0], n == 1, Which[ x < -1, 0, x < 0, x + 2*x^2 + x^3, x < 1, x - 2*x^2 + x^3, True, 0]]
myarg[ a_ ] := Module[ {rcd = a}, (* ローカル変数を定義・初期化 *) While[rcd >= Pi, rcd = rcd - 2 Pi]; While[rcd < - Pi, rcd = rcd + 2 Pi]; rcd];
Expand # 展開
Factor # 因数分解
Cancel # 約分
Simplify[ hoge, [assump]] # 式hogeを簡単にする。assumpで仮定を指定できる。 FullSimplify[ hoge, [assump]] # 初等関数と特殊関数を含む式hogeを簡単にする。
Expand[expr] 式の積とベキ乗の項を展開 ExpandAll[expr] すべての項にExpandを適用 Factor[expr] 因子の積の形に変換 Together[expr] 通分し単一分数にまとめる Apart[expr] 簡単な形の分母を持った複数の分数項に展開 Cancel[expr] 分母と分子の共通因子を約分する Simplify[expr] 式 exprに代数変形を試し,最も短い式の形を探す {Numerator[expr], Denominator[expr]} 式 exprの分子と分母のリストを作る Collect[expr,x] xの同じ次数の項をまとめる FactorTerms[expr,x] xに依存しない因子を取り出す TrigExpand[expr] 三角関数の式を項の和の形に展開 TrigFactor[expr] 三角関数の式を項の積の形に分解 TrigReduce[expr] 整数倍の角度を用いて三角関数を簡約 TrigToExp[expr] 三角関数を指数関数に変形 ExpToTrig[expr] 指数関数を三角関数に変形 FunctionExpand[expr] 特殊関数やその他の関数を展開 ComplexExpand[expr] すべての変数が実数からなるとした上で式を展開 PowerExpand[expr] 例えば(xy)^p を x^p y^p に展開する
(例) NN = 10; For[j = 0, j < NN, j++, a = 1 - ((j - 1/2*(NN - 1))/(1/2*(NN + 1)))^2; Print["j=", j, " a=", N[a, 3]]; ];