NewtonianMethod

Size: px

Start display at page:

Download "NewtonianMethod"

杻羿方
5 years ago
Views:

1 1 固有値問題入門 ( ニュートン法入門 ) 微分可能な方程式 = 0 の解を数値的に求めることを考えようその目的のためにニュートン法を説明しようこのときには複素正則であっても以下の議論は成り立つ今 = 0 の近似解として = が与えられているとしようそこで + = 0 となる補正を求めることを考えようそのために次のように + をテイラー展開する + = ここでの一次までの近似を採用するすると上式は + = 0 2 となるこれから = 3 となるそこで = + とおくと = 4 となるそこでをで置きかえて, 式 4 を遂次繰り返すことにより正確な解を求めようというのがニュートン法であるこのときに最初の値のの選択が重要である 3 式の関係を図示しておこう o = / = + =

2 例題 1. = 2 = 0 の解をニュートン法で求めようここで微係数は次の式より数値的に求めた方が一般に導関数が計算できないので一般的である 2 = 5 は h の精度で成り立つことに注意しておこうこの例題を解くために JAVA のプログラム SquareRootMain.java を提示しよう //======== SquareRootMain.java ======== class SquareRootMain{ // public static void mainstring args[]{ double x0=10;// -10; // 最初の値 ( 試行値 ) double h=1.0e-3*x0; double x=x0; double xf; double dx; double x1,x2; double fn,dfn; forint i=0;i<20;i++{ xf=x; x1=x+h; x2=x-h; fn=fnx; dfn=fnx1-fnx2/2.0*h;// 微係数 dx=-fn/dfn; h=1.0e-3*dx; x+=dx; // 漸化式 ifmath.absxf-x/xf<1.0e-11break; // 精度のための条件 System.out.println"i="+i+" x="+x; // static double fndouble x{ double y=x*x-4.0; return y; // //========================

3 3 例題微分方程式で定義される簡単な固有値問題を考えよう関数 yx は次の微分方程式で定義される = y 6 但し = 0= 0= ω を満足するものとするこのときに = を満足する ω の値を求めよよく知られているようにこの微分方程式の一般解は = + である = 0= より求める解は = である 1 = 0 より =0 であるから = = 1,2,3, である我々はこの問題の微分方程式を数値的に解きニュートン法で固有値を求め数値的に解決しようそしてこの学習から解析的に解けない量子力学の固有値問題を解く手立てが得られるそのために微分方程式微分方程式の数値解法についての数値解法についてまずまず述べておこう微分方程式 = 型の数値解法 =a で次のようにテイラー展開をすると + h + h 2 = h + 2 4! h + 同じように + h + h 2 = h + 2 4! h + 上式 2 4! h 下式で 2 4! h を消去 1 h h 10h + h + h + 1 h h

4 4 =! h + h の 5 次の精度で求めると 1 + h + h + 1 h h 2 + =0 (7) と h を与えると + h を 5 次の精度で次々と求めることができる動径部分のシュレーディンガー方程式 ħ +1ħ2 +[E Vr ]R=0 2 2 ところでは粒子の質量で,V r はポテンシャルである R= とおくと =r より数値計算で解れるシュレーディンガー方程式 ħ +[ ħ +[E Vr +1ħ E= ħ ] =0 2 E Vr ]=0 これは正に = 型である従って前述の方法で解くことが出来るそれでは具体的に例題 2. を解いて行こうまず微分方程式の解 y はとに依存するので解を, と書こうその時は = である解くべき方程式は, = であるそこで 7 の漸化式を使って [ ] の区間で微分方程式を数値的に解いて行こうそこで [ ] の区間を N 等分しその刻みを h とすると h=

5 5 そこで座標に関する配列を次のように定義するここで [0] = 0.0, [1] = h,, [] = 1.0, [ + 1] = h そしてに対応する解に関する配列を次のように定義する [0] = [0], [1] = [1],, [] = [], [ + 1] = [ + 1] そして 7 の漸化式の係数に対応する配列 c1 と c2 を次のように定義する 1[] = 1.0 [] 2[] = h2 12 [] すると 7 の漸化式は 1[ + 2]. [ + 2] + 1[]. [] 2[ + 1]. [ + 1] = 0 と書けてこれより [ + 2] = 2[ + 1]. [ + 1] 1[]. []/1[ + 2] 8 を遂次計算すれば解を求めることができるしかしこの際に初めの 2 項の値が必要である, をでテイラー展開すると初期条件から =0 +r 0+ = が 0 に近いところで近似すれば [0] =. [0], [1] =. [1] とすれば良いここでと, 1 の関係を表示する SinShapeMain.Java プログラムを提示しよう //================= SinShapeMain.Java =============== class SinShapeMain{ public static void mainstring args[]{ double dω=math.pi/20; double fn; System.out.println"a="+SolDifEq.a; SolDifEq ds0=null;

6 6 double ω=0.0; whileω<=math.pi*6.0{ ds0=new SolDifEqω; fn=ds0.sol0.0,soldifeq.a;// 0.0~1.0 System.out.printlnω+","+fn;//ω と yω,1 の関係を表示 ω+=dω; // class SolDifEq{ double ω; static double a=1.0; int N=1000;///////// 区間の刻み数 ///////////////////// double h; double[] xc=new double[n+2];//x 座標の配列 double[] c1=new double[n+2];// 漸化式の係数の配列 double[] c2=new double[n+2];// 漸化式の係数の配列 double[] ys=new double[n+2];// 解 y の配列 // SolDifEqdouble ω0{ this.ω=ω0; // double Soldouble x1,double x2{ h=x2-x1/n; double h2=h*h/12.0; double h10=h2*10.0; xc[0]=x1; forint i=0;i<n+1;i++{ xc[i+1]=xc[i]+h; c1[0]=1.0-h2*gxc[0]; c2[0]=2.0+h10*gxc[0]; forint i=0;i<n+2;i++{ c1[i]=1.0-h2*gxc[i];

7 7 c2[i]=2.0+h10*gxc[i]; ys[0]=ω*xc[0]; ys[1]=ω*xc[1]; forint i=0;i<n;i++{ ys[i+2]=c2[i+1]*ys[i+1]-c1[i]*ys[i]/c1[i+2];// 解 y の漸化式 double Sol=ys[N]; return Sol; // double gdouble r{ double y=-ω*ω; return y; // 計算結果よりと, 1 の関係のグラフを提示しよう 1.5 y(ω,1 ) ω y(ω,1 ) この図からの取り方により,2,3 に収束する例えば =7.0 に取れば 2 に収束するつぎに固有値とその時の微分方程式の解, 1 を求めるプログラム SinFn.java を提示しよう //============ SinFn.java======================== class SinFn{ public static void mainstring args[]{ double ω0=6.0;//;//+6.0;

8 8 double ω=ω0; double ωf; double h= *ω; double dω; System.out.println"a="+SolDifEq.a; SolDifEq ds1; SolDifEq ds2; SolDifEq ds0=null; forint i=0;i<40;i++{// yω,1=0 となる ω をニュートン法で求める double ω1=ω+h; double ω2=ω-h; ωf=ω; ds1=new SolDifEqω1; ds2=new SolDifEqω2; ds0=new SolDifEqω; double dfn=ds1.sol0.0001,soldifeq.a-ds2.sol0.0001,soldifeq.a/2.0/h;// 微係数 double fn=ds0.sol0.0001,soldifeq.a; System.out.println"i="+i+","+ω+","+fn; dω=-fn/dfn; h=0.001*dω; ω+=dω; ifmath.absωf-ω/ωf<1.0e-7break;// 精度 System.out.println" n 固有値 ω="+ω+" nr,yω,r"; forint i=0;i<ds0.n+1;i++{ ifi % 10 ==0{System.out.printlnds0.xc[i]+","+ds0.ys[i]; // class SolDifEq{// 微分方程式を解く double ω; static double a=1.0; int N=1000;////////// 刻み数 //////////////////// double h;// 刻み double[] xc=new double[n+2];

9 9 double[] c1=new double[n+2]; double[] c2=new double[n+2]; double[] ys=new double[n+2]; // SolDifEqdouble ω0{ this.ω=ω0; // double Soldouble x1,double x2{// x1 は始点座標 x2 終点座標 h=x2-x1/n; double h2=h*h/12.0; double h10=h2*10.0; xc[0]=x1; forint i=0;i<n+1;i++{ xc[i+1]=xc[i]+h; c1[0]=1.0-h2*gxc[0]; c2[0]=2.0+h10*gxc[0]; forint i=0;i<n+2;i++{ c1[i]=1.0-h2*gxc[i]; c2[i]=2.0+h10*gxc[i]; ys[0]=ω*xc[0]; ys[1]=ω*xc[1]; forint i=0;i<n;i++{ ys[i+2]=c2[i+1]*ys[i+1]-c1[i]*ys[i]/c1[i+2]; double Sol=ys[N];//yω,1 を返す return Sol; // double gdouble r{// gr double y=-ω*ω; return y; // //=================================== プログラム実行結果の一部を選択しメモ帳を開いて張り付けて txt ファイルに保存して

10 Excel でそのテキストファイルを開き挿入散布図でグラフ化すれば視覚化は容易である量子力学における固有値問題原子核内の核子の運動を 1 体ポテンシャルで記述するモデルがある独立粒子モデルと言われるものであるそこでの核子の従う波動方程式は次のように与えらる動径部分のシュレーディンガー方程式 ħ +[E Vr +1ħ ]R=0 9 ところでは粒子の質量で,V r は 1 体ポテンシャルである R= とおくと =r よりシュレーディンガー方程式は +[m 2 E Vr ħ ]=0 10 我々は 9 式を次の二つの境界条件のもとに E の固有値を求めたい 10 では r では 9 式よりは +1 2 =0 に支配されるここで = とおいて上式に代入すると = =0 より =, + 1 が得られる従って r= で正則な方の解を取らなくてはいけないので = である漸化式 (8) を解くときに必要な最初の 2 項はこの式 = を使うもう 1 度次の二つの境界条件をまとめると 1 r では = 11 2 では 12

11 11 もちろんこの微分方程式 10 は解析的には解けないので数値的に解く必要がある例題 3. A. ボーア B.R モッテルソンの共著原子核構造 Ⅰ 単一粒子運動 ( 講談社 ) 昭和 54 年の 195 ページを参照すると原子核内の核子の運動の1 体ポテンシャル Vrは次のように記述される = 0 +. 但し = [1 + ] fr 陽子数 Z, 中性子数 N 質量数 A とすれば 0= Mev = Mev = _0 / = 1.27 = 0.67 このときに +[ ħ 2 E Vr ] =0 11 から = 波のときのエネルギー固有値 E を求めようただし = 940, ħ = である考慮すべき境界条件は11と12よりと r では = 12 では 13 の 2 つである

12 我々の目的のためにはそのために微分方程式 (11) を数値的に充分に大きな r=r まで解いたときの (11) 式の解を, とすると r=r で 12, = 14 となる方程式を解いて E を決定すれば良いそのためにニュートン法を利用するここでと, の関係を求める java プログラム WoodsSaxonTrial.java を提示しよう //=========== WoodsSaxonTrial.java ================ class WoodsSaxonTrial{ public static void mainstring args[]{ ToSolve.A=208;///// 質量数 ////// ToSolve.Zn=82;///// 陽子数 ////// ToSolve.Nn=ToSolve.A-ToSolve.Zn;////// 中性子数 ToSolve.R=ToSolve.r0*Math.powToSolve.A,1.0/3.0;//// ポテンシャルの半径 R// System.out.println"A="+ToSolve.A+" Z="+ToSolve.Zn+" N="+ToSolve.Nn; double E0=-60.0;// エネルギーの下限 double l=0.0;//0.0;//// 粒子の軌道角運動量 /////// double j=0.5;//0.5;//// 粒子の全角運動量 /////// double E=E0; double h= *e; double de=1.0; System.out.println"R="+ToSolve.R+"fm"; ToSolve ds0=null; System.out.println" nemev,ye R"; whilee<0{ ds0=new ToSolveE,j,l; double fn=ds0.sol0.0001,4.0*tosolve.r;//x の [0.0001,4.0*ToSolve.R] の区間で解く System.out.printlnE+","+fn; E+=dE;

13 13 // class ToSolve{ double E; double j;//// 粒子の全角運動量 /////// double l;//// 粒子の軌道角運動量 /////// double m=940;// 核子の質量 in Mev double hc=197.3;//mev.fm static double a=0.67;//fm static double r0=1.27;//fm static double A; static double Zn; static double Nn; static double R; int N=1000;////////////////////////////// double h; double V0= *Nn-Zn/A; double Vls= *Nn-Zn/A; double[] xc=new double[n+2]; double[] c1=new double[n+2]; double[] c2=new double[n+2]; double[] ys=new double[n+2];// 微分方程式の解の配列 // ToSolvedouble E0,double j,double l{ this.e=e0; this.j=j; this.l=l; R=r0*Math.powA,1.0/3.0; //System.out.println"V0="+V0+" vls="+vls; // double Soldouble x1,double x2{ h=x2-x1/n; double h2=h*h/12.0; double h10=h2*10.0; xc[0]=x1; forint i=0;i<n+1;i++{ xc[i+1]=xc[i]+h;

14 14 c1[0]=1.0-h2*gxc[0]; c2[0]=2.0+h10*gxc[0]; forint i=0;i<n+2;i++{ c1[i]=1.0-h2*gxc[i]; c2[i]=2.0+h10*gxc[i]; ys[0]=xc[0]; ys[1]=xc[1]; forint i=0;i<n;i++{ ys[i+2]=c2[i+1]*ys[i+1]-c1[i]*ys[i]/c1[i+2]; double Sol=ys[N]; return Sol; // double gdouble r{ double y=-2.0*m/hc*hc*e-vr+l*l+1.0/r*r; return y; // double Vdouble r{ // ポテンシャル Vr double y=v0*fr+vls*r0*r0*dfr/r*0.5*j*j+1-l*l ; return y; // double fdouble r{ //fr double y=1.0/1.0+math.expr-r/a; return y; // double dfdouble r{ //frの微分 double yn=math.exp1.0; double y=-1.0/a*yn*math.pow1.0+yn,-2.0; return y; //============================================= ここでと, の関係のグラフを示しておこうここでは ye,rが, に対応している

15 15 log(abs(y(e,r))) s 2s 3s 4s E0=-40Mev E0=-30Mev E0=-15Mev E0=-2Mev E0 はニュートン法における最初の試行値である 4 2 E (Mev) log(abs(y(e,r))) つぎにこの例題の固有値とその時の固有関数を計算する JAVA プログラム WoodsSaxon.java を提示しよう //============== WoodsSaxon.java =================== class WoodsSaxon{ public static void mainstring args[]{ DifSolution.A=208;/////////// DifSolution.Zn=82;/////////// DifSolution.Nn=DifSolution.A-DifSolution.Zn; DifSolution.R=DifSolution.r0*Math.powDifSolution.A,1.0/3.0; System.out.println"A="+DifSolution.A+" Z="+DifSolution.Zn+" N="+DifSolution.Nn; double E0=-40.0;/////// 試行値 E01s=-40.0* E02s=-40.0* E03s=-15.0* E04s=-2.0 double l=0.0;/////////// double j=0.5;/////////// double E=E0; double h= *e; double de; System.out.println"R="+DifSolution.R; DifSolution ds1; DifSolution ds2; DifSolution ds0=null; System.out.println"i 回数,E エネルギー固有値 ";

16 16 forint i=0;i<40;i++{ double E1=E+h; double E2=E-h; ds1=new DifSolutionE1,j,l; ds2=new DifSolutionE2,j,l; ds0=new DifSolutionE,j,l; double dfn=ds1.sol0.0001,3.0*difsolution.r-ds2.sol0.0001,3.0*difsolution.r/2.0/h; double fn=ds0.sol0.0001,3.0*difsolution.r; System.out.println"i="+i+","+E+","+fn; ifdouble.isnan dfn dfn==0.0break; de=-fn/dfn; h=0.001*de; E+=dE; double norm=ds0.norm_shimpsonds0.ys; forint i=0;i<ds0.n+1;i++{ds0.nys[i]=norm*ds0.ys[i]; System.out.println"Norm="+ds0.norm_shimpsonds0.ys+" NNorm="+ds0.norm_shimpsonds0.nys; System.out.println"x, 計算波動関数 ys, 規格化波動関数 ysnormed_ys"; forint i=0;i<ds0.n+1;i++{ ifi % 10 ==0{System.out.printlnds0.xc[i]+","+ds0.ys[i]+","+ds0.nys[i]; // class DifSolution{ double E; double j; double l; double m=940; double hc=197.3; static double a=0.67; static double r0=1.27;

17 17 static double A; static double Zn; static double Nn; static double R; int N=1000;/////////// 区間の刻み数 /////////////////// double h;/////////// 区間の刻み /////////////////// double V0= *Nn-Zn/A; double Vls= *Nn-Zn/A; double[] xc=new double[n+2]; double[] c1=new double[n+2]; double[] c2=new double[n+2]; double[] ys=new double[n+2]; double[] nys=new double[n+2];// 規格化波動関数の配列 // DifSolutiondouble E0,double j,double l{ this.e=e0; this.j=j;// 全角運動量 this.l=l;// 軌道角運動量 R=r0*Math.powA,1.0/3.0; //System.out.println"V0="+V0+" vls="+vls; // double Soldouble x1,double x2{ //double x1,double x2 は始座標終座標 h=x2-x1/n; double h2=h*h/12.0; double h10=h2*10.0; xc[0]=x1; forint i=0;i<n+1;i++{ xc[i+1]=xc[i]+h; c1[0]=1.0-h2*gxc[0]; c2[0]=2.0+h10*gxc[0]; forint i=0;i<n+2;i++{ c1[i]=1.0-h2*gxc[i]; c2[i]=2.0+h10*gxc[i];

18 18 ys[0]=xc[0]-math.powxc[0],3; ys[1]=xc[1]-math.powxc[1],3; forint i=0;i<n;i++{ ys[i+2]=c2[i+1]*ys[i+1]-c1[i]*ys[i]/c1[i+2]; double Sol=ys[N]; return Sol; // double gdouble r{ double y=-2.0*m/hc*hc*e-vr+l*l+1.0/r*r; return y; // double Vdouble r{ double y=v0*fr+vls*r0*r0*dfr/r*0.5*j*j+1-l*l ; return y; // double fdouble r{ double y=1.0/1.0+math.expr-r/a; return y; // double dfdouble r{ double yn=math.exp1.0; double y=-1.0/a*yn*math.pow1.0+yn,-2.0; return y; // double norm_shimpsondouble[] xys{ // シンプソンの積分公式 int num=xys.length; int xnum=num-2/2; double s=0.0; forint i=0;i<xnum;i++{ s+=h/3.0*xys[2*i]*xys[2*i]+4.0*xys[2*i+1]*xys[2*i+1]+xys[2*i+2]*xys[2*i+2]; return 1.0/Math.sqrts; //

19 19 //=========================================== WoodSaxon.java により計算された結果を提示しておこう WoodSaxon.java により計算された結果状態試行値 E0(Mev) 始点 x1(fm) R(fm) 終点 x2(fm) 繰返回数収束値 (Mev) ( 計算固有値 ) 1 / (=3R) / (=3R) / (=3R) / (=3R) WoodSaxon.java により計算された規格化固有関数 dr=1 6.00E-01 1S 5.00E-01 2S 4S 4.00E E-01 3S 2.00E E E E E E E E E E E E E E E-01 1s 2s 3s 4s 例題 4. 一次元調和振動子に関する固有値問題一次元調和振動子のシュレーディンガー方程式は次のようになる ħ K = 15 ここで =/ とおき = m ħ + = 0 とおいて = m x とおくと ħ 16 となる但し = である ħ そこで境界条件 ± = もとに固有値もとに固有値を求めよ

20 20 ここで 16 式は ± = 0 の境界条件のもとでは解は一意に決まることを示しておこうそこで 16 式に 2 つの解との二つがあったとしてみようすると 16 より = 0 17 = 0 18 の二つが成り立つそこでを作ると = d d =0 を得る = 定数 19 が成り立つこれより従って 19 は定数 = ± ± ± ± = 0 =0 20 ここで = と置いて = を 20 に代入すると + = =0 なるので 0 で =0 となる従って 0 なる適当な開区間で C= 定数であるここで C 0 となることを示そうもし C=0 となれば 0 なる適当な開区間で恒等的に 0 これより =0 かつ = 0 より 16 の解が自明的な解 0 ではないとすれば C 0 となるつぎに示されている解の一意性を使うことにより = の関係を 0 なる適当な開区間から定義域全体に拡張することによって 16 の固有値は縮退していないという重要な結果が導かれるそのために解の一意性を示しておく = 0, かつ = 0 のとき微分方程式 16 の解は必ず 0 となることを証明しよう ξ を ξ < ξ になるようにしかし任意にととり閉区間 ξ ϵ [ξ, ξ ] おけるの最大値を K m= + とおき m の最大値を M とする

21 21 17 即ち, m= + である従って = = = より d d d + d d m 1 + m d 19 d d d + d 18 m 1 + m これを (19) の右辺に代入すると d 1 + m 1 + d =1 +! m 1 + 2! またこれを19の右辺に代入しこの操作を繰り返すとをえる m [ ]! 20 ここで n とすると m=0 なるから 0 かつ 0 となる ξ > ξ 同じように示すことが出来るこのことから定義域内の点で =, = が成り立つとき定義域内の点で恒等的にが成り立つ ( 解の一意性 ) 16 の解の一つを u とすれば 16 の解はすべて u に一次従属である即ち 16 の任意の解をとすれば =C 定数 u と書ける量子力学的に言えば 16 は即ち 15 は縮退をしていないない 16おいて = で

22 22 よりをで書き直すと + = 0 + = 0 からも 16 の解であるから定数を使って解の従属性より = 21 と書けるをで置きかえると = が得られ = = より =1 が得られる従って 21 から = ± となり 16 の解は偶関数か奇関数に限られることがわかる我々の目的のためには微分方程式 (16) を =0から数値的に充分に大きな = まで解いたときの (16) 式の解を, とするとき = で, = 0 22 となる方程式を解いてを決定すれば良いそのためにニュートン法を利用するこのときに 7 の関係の漸化式を使い解いて行くわけであるがそのためには漸化式の最初の 2 項が必要であるここで 16 の解は偶関数と奇関数に限られることを考える 7 の刻みを h とすると解を原点 =0 のまわりでテイラー展開すると h=0 + h 0 + h が奇関数のときはより h=h

23 23 0 = 0 h = 任意の定数 23 と最初の 2 項を決めればよいここでの任意の定数倍もまた解である事実を用いたつぎにが偶関数のときはより h=0 + h = 0 h 2 h 0 + 0= 任意の定数 h= と最初の 2 項を決めればよいここでもの任意の定数倍もまた解である事実を用いたところで < 0 のときのの値は 0 のとき漸化式で求めた値用いて偶関数あるいは奇関数の性質を使って拡張する JAVA による一次元調和振動子に関する固有値問題のプログラムは偶関数用には HarmonicOscillaterEven.java に奇関数用には HarmonicOscillaterOdd.java に記載されているが 16 の固有値の分布がどうなるかを知るために λ と, の関係のグラフを固有関数が偶関数の時と奇関数の時の 2 つを示しておこうこのグラフから固有値の試行値を初めにどう取ればよいかがわかるこれら 2 つのグラフは HarmonicOscillaterEvenShape.java と armonicoscillateroddshape.java の計算結果より Excel を使って作図されたグラフより偶関数のときは 1,5,9,13, 奇関数のときは 3,7,11,15, が固有値になっていることがわかる

24 24 Log(Abs(u(λ X))) 18 偶関数のとき λ Log(Abs(u(λ,X))) 20 奇関数のとき λ

目次第 1 章はじめに 1 1 本研究の目的 1 1 2 先行研究 2 1 3 本研究の構成 9 1 4 略語収集の方針 10 第 2 章原形分割単位 (1): 段 2 1 略語構成過程における原形分割 13 2 2 段の

目次第 1 章はじめに 1 1 本研究の目的 1 1 2 先行研究 2 1 3 本研究の構成 9 1 4 略語収集の方針 10 第 2 章原形分割単位 (1): 段 2 1 略語構成過程における原形分割 13 2 2 段の日本大学大学院文学研究科学位論文現代中国語略語研究構成過程上の制約と数詞の位置決定規則学位申請者星健一 Kenichi HOSHI 平成 25 年 11 月 7 日目次第 1 章はじめに 1 1 本研究の目的 1 1 2 先行研究 2 1 3 本研究の構成 9 1