計算物理学２

第２回演習

ver. 2018/4/19
以下の例題を Fortran もしくは C 言語を用いてプログラム化せよ。
四則演算(1) 標準入力から実数 a と b を読み込み、その和、差、積、商を標準出力に出力。
・b にゼロを代入するとどうなりますか。
DO ループ

(2)

NX

k=1
<latexit sha1_base64="RbrC/O/+NaLXA/c8zqnfpKtThIQ=">AAAClXichVHLSsNAFD3GV62PRl0ouCkWxVW5EUERhKIirrQ+WoWqJYmjhuZFkha0+AP+gAvdKFQQP8APcOMPuPATxKWCGxfepgFRUe8wM2fO3HNnDldzTcMPiB6bpOaW1rb2WEe8s6u7JyH39uV9p+zpIqc7puNtaqovTMMWucAITLHpekK1NFNsaKW5+v1GRXi+4djrwaErti113zb2DF0NmCrK8pZftorV0oxyvLOULO2MF+UUpSmM5E+gRCCFKLKOfIst7MKBjjIsCNgIGJtQ4fMoQAHBZW4bVeY8RkZ4L3CMOGvLnCU4Q2W2xOs+nwoRa/O5XtMP1Tq/YvL0WJnECD3QNb3QPd3QE73/Wqsa1qj/5ZB3raEVbjFxMrj29q/K4j3AwafqD4XG2X97CrCHqdCLwd7ckKm71Bv1K0enL2vTqyPVUbqkZ/Z3QY90xw7tyqteWxGrZ4hzg5Tv7fgJ8uNphdLKykQqMxu1KoYhDGOM+zGJDBaRRY7freAcNVxJA9KMNC8tNFKlpkjTjy8hLX8AtsaV5g==</latexit>

k2

NXおよび k=1
<latexit sha1_base64="YMXOWn6xmmnzzqFNBiIE1FqLLsI=">AAAClXichVHLSsNAFD2N7/qKulBwUyyKq3KjgiIUioq40vqoCj5KEkcNzYskLWjxB/wBF7pRqCB+gB/gxh9w0U8QlwpuXHibBkRFvcPMnDlzz505XM01DT8gqsakhsam5pbWtnh7R2dXt9zTu+47RU8XOd0xHW9TU31hGrbIBUZgik3XE6qlmWJDK8zW7jdKwvMNx14LjlyxY6kHtrFv6GrAVF6Wt/2ilS8X0srJ7mKisDuel5OUojASP4ESgSSiyDryHbaxBwc6irAgYCNgbEKFz2MLCgguczsoM+cxMsJ7gRPEWVvkLMEZKrMFXg/4tBWxNp9rNf1QrfMrJk+PlQkM0yPd0As90C090fuvtcphjdpfjnjX6lrh5rtPB1bf/lVZvAc4/FT9odA4+29PAfYxFXox2JsbMjWXer1+6fjsZXV6Zbg8Qlf0zP4uqUr37NAuveqVZbFyjjg3SPnejp9gfSylUEpZnkhmZqJWtWIQQxjlfkwigwVkkeN3S7hABddSv5SW5qT5eqoUizR9+BLS0ge45pXn</latexit>

k3の N=1 から 30 までの値。

(3) N! (階乗)の値(N が 1 から 30 まで)
・N がいくつまで正確に計算できますか。
(4) 整数 N=1,...30 に対して単精度と倍精度で N の値を出力(sqrt(x)関数を使う)。
配列(5) 行列 A(1:3,1:3)と B(1:3,1:3)の積の計算(Fortran の場合は MATMUL を使う場合と使わない場合二通り)。
ただし A の行列要素 aij は aij = e<latexit sha1_base64="vUrgH5SzVTHWw9hQmPxsic5JTWE=">AAACmHichVHLLgRBFD3ae7wGG2EzMSE2JrdFQiSSCQvsGAaJx6S71VD0dHe6eyYZbX7AD1hYCImE+AAfYOMHLHyCWJLYWLjT04kguJWqOnXqnlt1cnXHlJ5P9Fin1Dc0NjW3tMba2js6u+LdPaueXXQNkTVs03bXdc0TprRE1pe+KdYdV2gF3RRr+sFs9X6tJFxP2taKX3bEVkHbtWReGprPVC7eq+UCuV9JTCfEdjB6JEf3jyq5eJJSFEbiJ1AjkEQUi3b8FpvYgQ0DRRQgYMFnbEKDx2MDKggOc1sImHMZyfBeoIIYa4ucJThDY/aA110+bUSsxedqTS9UG/yKydNlZQJD9EDX9EL3dENP9P5rrSCsUf1LmXe9phVOruu4b/ntX1WBdx97n6o/FDpn/+3JRx6ToRfJ3pyQqbo0avVLhycvy1OZoWCYLuiZ/Z3TI92xQ6v0alwuicwpYtwg9Xs7foLVsZRKKXVpPJmeiVrVggEMYoT7MYE05rGILL9bxhmucK30K2llTlmopSp1kaYXX0LJfADfS5dE</latexit>|i j|bij は bij = e<latexit sha1_base64="PiOqhrUvI21cBeQNC1ugrQdDPTI=">AAACl3ichVHLLgRBFD3ae7wGG2IzMSE2JrdFQiRCSMTSYJB4TLpbDUVPd6e7Z4I2P+AHLCQSEkF8gA+w8QMWPkEsSWws3OnpRBDcSlWdOnXPrTq5umNKzyd6rFKqa2rr6hsaY03NLa1t8faOJc8uuIbIGLZpuyu65glTWiLjS98UK44rtLxuimV9d7p8v1wUridta9Hfd8R6XtuyZE4ams9UNt6hZwO5U0qMJ8RGcCgHdw5L2XiSUhRG4idQI5BEFHN2/BZr2IQNAwXkIWDBZ2xCg8djFSoIDnPrCJhzGcnwXqCEGGsLnCU4Q2N2l9ctPq1GrMXnck0vVBv8isnTZWUCffRA1/RC93RDT/T+a60grFH+yz7vekUrnGzbUdfC27+qPO8+tj9Vfyh0zv7bk48cRkMvkr05IVN2aVTqFw+OXxbG5vuCfjqnZ/Z3Ro90xw6t4qtxkRbzJ4hxg9Tv7fgJloZSKqXU9HBycipqVQN60IsB7scIJjGLOWT43T2c4hJXSrcyocwos5VUpSrSdOJLKOkPVseXDg==</latexit>
sha1_base64="dO1RVWgbv0QBXFnWGIM/DJuxBX8=">AAACdXichVG9TsJQGD1UVEQU3ExciARjHMhXHTROJi6O/MhPgoS05YINpW3aQoLEF2B1cHDSxMH4AD6Aiy/gwCMYR0xcHPwoJEaJ+jXtPffc73y3J0e1Dd31iAYBaSY4OzcfWggvRsJLy9FYpOBabUcTec0yLKekKq4wdFPkPd0zRMl2hNJSDVFUm4ej82JHOK5umcde1xaVltIw9bquKR5T6WosQSnyKz4N5AlIYFJW7AEnqMGChjZaEDDhMTagwOWnDBkEm7kKesw5jHT/XOAcYda2uUtwh8Jsk78N3pUnrMn70UzXV2t8i8Gvw8o4kvRMdzSkJ7qnF/r4dVbPnzH6ly6v6lgr7Gq0v5p7/1fV4tXD6ZfqD4XK3X978lDHnu9FZ2+2z4xcauP5nbPLYW4/m+xt0A29sr9rGtAjOzQ7b9ptRmSvEOZ85J9pTIPCdkqmlJwhhLCGdWxyDLs4wBHSyPN1NfRxIQWlLUke5ygFJoGu4FtJO59jaI2O</latexit>
sha1_base64="EAg32O9wieweTtupxTmfNBQYJzM=">AAACjHichVHNLgNRGD3Gf/2VFbFpNBUbzTc2RCIkEumyRZFUNTPj4jKdmcxMG0z7Al7AQiIhEcQDeAAbL2DhEcSSxMbC12kTQfBNZu65537nu3NydMeUnk/02KQ0t7S2tXd0Rrq6e3r7ov3dK55dcg2RNWzTdtd0zROmtETWl74p1hxXaEXdFKv63nztfLUsXE/a1rJ/4Ih8Udu25JY0NJ+pQnRALwRytxqbiYmNoCLHdyvVQjROSQor9hOoDRBHo9J29Bbr2IQNAyUUIWDBZ2xCg8dPDioIDnN5BMy5jGR4LlBFhLUl7hLcoTG7x99t3uUarMX72kwvVBt8i8mvy8oYEvRA1/RC93RDT/T+66wgnFH7lwNe9bpWOIW+o8Glt39VRV597Hyq/lDo3P23Jx9bmAq9SPbmhEzNpVGfXz48flmaXkwEo3ROz+zvjB7pjh1a5VfjIiMWTxDhgNTvcfwEKxNJlZJqhtCBYYxgjGOYxBxSSCPL1+3jFJe4UoaUWWWhHqXS1Mh0AF9KSX0AzN2WCw==</latexit>
sha1_base64="oEZ5gg4a1MpEvUsaLF+dhsjragE=">AAACl3ichVHLLgRBFD3aa4zXYENsJibExuS2DZEIIRFLg0HiMeluhTI93Z3ungna/IAfsJBISATxAT7Axg9Y+ASxJLGxcKenE0FwK1V16tQ9t+rk6o4pPZ/osUapratvaIw1xZtbWtvaEx2dS55ddA2RNWzTdld0zROmtETWl74pVhxXaAXdFMt6frpyv1wSridta9Hfd8R6Qdu25JY0NJ+pXKJTzwVyt5wcT4qN4FAO7R6Wc4kUpSmM5E+gRiCFKObsxC3WsAkbBoooQMCCz9iEBo/HKlQQHObWETDnMpLhvUAZcdYWOUtwhsZsntdtPq1GrMXnSk0vVBv8isnTZWUS/fRA1/RC93RDT/T+a60grFH5yz7velUrnFz7UffC27+qAu8+dj5Vfyh0zv7bk48tjIZeJHtzQqbi0qjWLx0cvyyMzfcHA3ROz+zvjB7pjh1apVfjIiPmTxDnBqnf2/ETLA2nVUqrGUpNTkWtiqEXfRjkfoxgErOYQ5bf3cMpLnGl9CgTyowyW01VaiJNF76EkvkAVYeXCg==</latexit>|i j|とする。

(6) n が 5 までのルジャンドル多項式の直交性を以下の方法で評価しなさい。
まず区間[-1,1]を N-1 等分し、それぞれの座標を要素数 N の配列 x(1:N)にセットする。端点を x(1)=-1, x(N)=1 として残りの点 x(2), x(3), ...., x(N)も求めて値を配列にセット。
配 列 P(1:N,0:5) を 用 意 し 、 同 じ [-1,1] 区 間 で ル ジ ャ ン ド ル 多 項 式 P0(x), P1(x),P2(x),P3(x),P4(x),P5(x)の値をセット。(例えば x=x(i)でのルジャンドル多項式 P2(x(i))の値を配列 P(i,2)に代入) ただし n=0,1,2,3,4,5 のルジャンドル多項式は以下で与えられる。

1

P0 (x) =1P1 (x) =x
1P2 (x) = (3x2 1)
21
P3 (x) = (5x3 3x)
21
P4 (x) = (35x4 30x2 + 3)
81
P5 (x) = (63x5 70x3 + 15x)8<latexit sha1_base64="wrH3Ff+mQ0YK39etlN1YW5IO/K8=">AAADMXichVHLTttQEB27tKVuKSndVGITNQIlQonmxkmJKlXisWEZHgEkApZ9e0MtHNuyncjUyg/0B7pgRaUuqnbXBR/QDT/AIruqu4olSGxYdPyQKh6lY9l35sw553o0hmuZfoA4lOR7I/cfPBx9pDx+MvZ0PPdsYt13eh4XLe5Yjrdp6L6wTFu0AjOwxKbrCb1rWGLD2FuM+xt94fmmY68F+67Y7uq7ttkxuR4QpOW+NTUshqX8m+k8a7eVpsayKkyqala1O57OIzaIqoOiGu5Uy6yU9NVb+vVwRy2rYUqoXSc0yIAYtbKK5DOjprT6LbRXdFG9PEs0dYbVw5KiaLkCVjCJ/M2EZUkBsmg6uSNow1twgEMPuiDAhoByC3Tw6dkCBgguYdsQEeZRZiZ9AQNQSNsjliCGTugefXep2spQm+rY00/UnG6x6PVImYcpPMEveIbH+BV/4+U/vaLEI/6XfTqNVCtcbfzDi9WL/6q6dAbw7q/qDoVB7LtnCqADjWQWk2ZzEySekqf+/fcfz1Zfr0xF0/gJT2m+QxziD5rQ7p/zz8ti5QDiBbHr67iZrFcrDCtsuVaYW8hWNQqT8BKKtI9ZmIMlaEILuDQhNaR5aUH+Lg/ln/KvlCpLmeY5XAn59A8m1sJf</latexit>

ルジャンドル多項式の配列を用いて、区分求積法で次の積分を計算する。
Z 1

Inm =1<latexit sha1_base64="6F4N5B0OfOV6VoA7bi1pwuFVRtI=">AAAC5HichVFRSxRRFP6csmzTXOsl6GVwMdYHlzMRFIIg1UO9xKqtCq4OM+NduzhzZ5i5u6wN+weit4gegqCgh+gH+ANE6A/04E8QH1V66aEzswOiop7L3PPd737fuXM4buTLRBPtDRjXrg/euDl0q3R7eOTOaHns7mIStmNPNLzQD+Nl10mEL5VoaKl9sRzFwglcXyy5m8+z+6WOiBMZqjd6KxKrgbOhZEt6jmbKLtuv7FQFPXPGbEql7XTK6q1ZZt1W1e4kpyBL612zmciAt3Zgp3KGJelrFpp9XVVOFtIcZdF8IXztmN1SyS5XqEZ5mOeBVYAKiqiH5W00sY4QHtoIIKCgGftwkPBagQVCxNwqUuZiRjK/F+ihxN42qwQrHGY3ed/g00rBKj5nNZPc7fErPn8xO01M0B/6SYf0m37RPv27sFaa18j+ZYuz2/eKyB59f3/h75WugLPG2xPXJQ6X1Zf3pNHC07wXyb1FOZN16fXrd959PlyYnp9IH9J3OuD+vtEe7XCHqnPk/ZgT81+QDcg6O47zYPFRzaKaNfe4MvusGNUQHmAcVZ7HE8ziJepo8Lu7OMARjo2W8cH4aHzqS42BwnMPp8L4+h9CDLA1</latexit>

Pn (x)Pm (x)dx ⇠NX1Pn (x(i))Pm (x(i)) xi=1

Δx は N から計算できる。Inm を 6x6 の行列として出力。
・N の値を 100, 1000,10000 と変えるとどうなるか。
・区分求積の和を N-1 までにした時と N までにした時で何が変わるか。
(7) 標準入力から２つの角度 a と b を実数で読み取り(単位：度)、cos(a+b)と sin(a+b)
を 3 通りの方法で計算。
1. cos(a+b), sin(a+b)をそのまま出力(ラジアンに変換)。
2. 角度 a 回転させる 2 次元の行列、角度 b 回転させる行列を配列を使って作り、それらの行列の積を計算して cos(a+b)と sin(a+b)に対応する成分を出力。ただし角✓cos ✓度θ回転させる回転行列は sin ✓<latexit sha1_base64="RTKqxPM5n6VVbTCs+/RZWd9frvY=">AAAC13ichVGxTtxAEH04kJADwhGaSGmsnDil4TRGSEFUKGko4cgBEUYn2yzHCntt2XsnyAmlQ4iOioIKpBRRPiD0NPwABQ1tQClBSpOCOZ9FSBAw1npn3r43u0/jRr5MNNFph/Gks+vps+7nuZ7evhf9+YGXs0lYjz1R8UI/jOddJxG+VKKipfbFfBQLJ3B9MeeufmidzzVEnMhQfdTrkVgMnJqSy9JzNEPVfNl2RU2qZhQ4OpZrG2bO9sLE1itCO2bRHLYTqbLKts1bVdG8RbSFWrppUc0XqERpmHcTK0sKyGIqzP+AjSWE8FBHAAEFzbkPBwl/C7BAiBhbRJOxmDOZngtsIMfaOrMEMxxGV/lf42ohQxXXrZ5Jqvb4Fp9XzEoTQ3RC3+iSjuk7XdCfe3s10x6tt6zz7ra1Iqr2b7+a+f2oKuBdY+Wv6gGFy+yHPWksYyz1ItlblCItl167f+Pz7uXMeHmoWaQD+sX+9umUjtihalx5X6dFeQ85HpD1/zjuJrMjJYtK1vRoYeJ9NqpuvMYbvOV5vMMEJjGFCt97iDP8xLnxyfhibBpbbarRkWkG8U8YO9fo3bFR</latexit>ia◆sin ✓cos ✓ で与えられる。

ib

3. 複素数 e と e の積を計算し、その実部と虚部から出力。

2