● 前回の講義のまとめ
★ ルンゲクッタ法の安定性
• 計算スキームを
x′(t) =λx(t), x(0) = 1 (1)
に適用したとき, λ <0 ならば, 数値解{xn} は xn →0 (n→ ∞)が成り立つべきである. これが 成り立つかどうかは, 計算スキームとz =λh に依存する. これを線形テスト問題と呼び,xn →0, n→ ∞が成り立つとき,その計算スキームはzで絶対安定であるという.
• 前進オイラー法を (1)に適用すると,xn= (1 +λh)n であるので,−2< z <0 の時に絶対安定とな る. λ=−1の時を考えれば,h≥2ととると, (1)の数値解は安定ではなくなる.これは,長時間の計 算を行う際にもhをあまり大きくとってはいけないことを意味している.
• 一般のp段陽的ルンゲクッタ法(p= 1, . . . ,4)の場合, xn+1 =Rp(z)xn
と書くことができるので,|Rp(z)|<1 となる領域
Rp={z∈C:|Rp(z)|<1}
となる領域（絶対安定領域）の内部に zが入れば絶対安定となる.
• 具体的にR(z)を計算すれば,
R1(z) = 1 +z, R2(z) = 1 +z+z2/2,
R3(z) = 1 +z+z2/2 +z3/3!, R4(z) = 1 +z+z2/2 +z3/3! +z4/4!
となり,任意のλ <0に対して,十分小さい h >0 を取れば絶対安定となる. また,ルンゲクッタ法
の場合には,段数が増加するほど絶対安定領域は広がる傾向がある.
★ シンプレクティック法
• 関数H: R2n −→R2n ((q, p)7→H(q, p))が与えられたとき, dq
dt =∇pH(q, p), dp
dt =−∇qH(q, p)
(2)
で定められる常微分方程式の初期値問題を考える.
Example 1 H(q, p) = 12(q2+p2)とおくと, dq dt =p, dp dt =−q
Example 2 H(q, p) = 12p2−cos(q)とおくと, dq dt =p, dp
dt =−sin(q)
となる. この方程式は,x′′(t) =−sin(x(t))を書き直したものである.
• 方程式(2) の一つの解を(q(t), p(t))と書くと, d
dtH(q(t), p(t)) = 0
が成り立つ. すなわち,H は解にそって不変な関数になる. このH をハミルトニアンと呼ぶ.
• このようなハミルトン系の方程式を陽的ルンゲクッタ法で解くと,ハミルトニアンの値が増大または 減少し,長時間にわたって誤差の少ない数値解を得ることが難しい.
• † 2n次正方行列M が
tM JM =J, J = 0 1
−1 0
!
をみたすとき,M はシンプレクティック行列であるといい,シンプレクティック行列全体は群をなし, それをSp(n)と書く. いま,M ∈Sp(n)ならばdet(M) = 1が成り立つ.
また,写像f: R2n −→R2nがシンプレクティックであるとは,任意の(q, p)∈R2nに対して,Df(q, p)∈
Sp(n)となることである.
• † 微分方程式(2) に対して,写像StH:R2n−→R2n を
StH(q0, p0) = (q(t), p(t)), (q(t), p(t))は,初期条件(q0, p0)に対する(2)の解
とおく. このとき,任意のtに対してStH はシンプレクティックとなる. すなわち,StH は面積保存写 像である.
• 方程式(2) の数値解法として,
qn+1=qn+h∇pH(qn, pn),
pn+1=pn−h∇qH(qn+1, pn) (3)
を考える. これをシンプレクティックオイラー法と呼ぶ. 一般のH に対しては,陰的解法となるが, H(q, p) =T(p) +U(q)
と書けているときには,
qn+1=qn+h∇pT(pn),
pn+1=pn−h∇qU(qn+1) (4)
と書け,陽的解法とすることができる.
• シンプレクティックオイラー法は,局所離散化誤差O(h2)となる. したがって1 次の公式である.
• † シンプレクティックオイラー法によって定まる写像SEh:R2n−→R2n を SEh(qn, pn) = (qn+1, pn+1)
とおくと,SEh はシンプレクティック写像となる. 実際, Example 1の方程式に対しては,線形写像と なり,
1 h
−h 1−h2
!
で定められる写像となる. この行列はSp(n)に属することは容易に確められる.
• † 写像SQh,SPh を
SQh: q p
!
7→ q+h∇pH(q, p) p
! ,
SPh: q p
!
7→ p
p−h∇pH(q, p)
! ,
と定めると,
SEh =SPhSQh が成り立つ.
DSQh = 1 h∇pH
0 1
!
, DSQh = 1 0
−h∇qH 1
! ,
であるので,SQh,SPh はともにシンプレクティック写像となる.
• † より高次の公式を得るには,
(SEp)h:=SPcpSQbp· · ·SPc1SQb1
を考える. この公式がp次となるように{bi}pi=1,{ci}pi=1 を定めることができる. 実際, p= 2, (St¨ormer-Verlet)
1 2
bi 0 1
ci 1/2 1/2
p= 3, (Ruth)
1 2 3
bi 7/24 3/4 -1/24
ci 2/3 -2/3 1
という公式を得ることができる.
• シンプレクティック法(SEn)h は,厳密な不変量
H(q, p) =˜ H(q, p) +hnHn(q, p) +· · · を持つことが知られている. 特に, Example 1の場合には,
H˜(q, p) =q2+p2+hqp が不変量となり,この量が定める楕円上に拘束される.
• したがって,ハミルトン系の方程式を解く際には,シンプレクティック法を用いることにより,長時間
• † 陰的ルンゲクッタ法は,パラメータ{aij},{bi},{ci}が
mij :=biaij+bjaij−bibj= 0, 1≤i, j ≤s をみたせば,シンプレクティックとなることが知られている.
• 1 段 の陰的ルンゲクッタ法では,
b1= 1, a11b1= 1 2
とパラメータを決めることができる. 実際,a11=c11= 1/2,b1= 1 とすればよい. したがって,この ルンゲクッタ法は2 次公式となり,シンプレクティックとなる. 具体的に書けば,
X=xn+f(tn+h/2, X), xn+1=xn+f(tn+h/2, X) となる. これを陰的中点法と呼ぶ.
• † 一般にp 段の陰的ルンゲクッタ法で, パラメータを適切に定めることにより, 2p次のシンプレク ティクなスキームを得ることができる.
● 講義資料
★
C
で行列を表現する方法1. double型の二重配列で表現する. （多分, もっとも標準的な方法）
【長所】 わかりやすい.
【短所】 関数に配列のサイズを明示的に指定する必要がある.関数に行列のサイズをわたさないと, 無駄な計算をすることになる.
2. double型の配列で表現する.
【長所】 関数定義が配列サイズと無関係になる.
【短所】 呼び出された関数側では,どこまでが１行かを知る方法がないため,関数に行列のサイズを 渡す必要がある.
3. 行列要素のデータと行列サイズを構造体にする.
【長所】 関数定義が単純になる.
【短所】 初期化が面倒？
実際の以下のプログラムは,紙面の関係上,空行などを取り除き,多少読みにくくなっている.
★ 実際のプログラム
☆ double 型の二重配列で表現する
/* $Id: matrix_0.c,v 1.2 2006-06-15 19:55:54+09 naito Exp $ */
#include <stdio.h>
#define MAX (10)
void print_matrix(double [][MAX], int, int) ;
void add_matrix(double [][MAX], double [][MAX], double [][MAX]) ; int main(int argc, char **argv)
{
double a[MAX][MAX] = {{1.0, 2.0}, {2.0, 3.0}} ; double b[MAX][MAX] = {{1.0, 2.0}, {2.0, 3.0}} ; double c[MAX][MAX] ;
print_matrix(a, 2, 2) ; print_matrix(b, 2, 2) ; add_matrix(c, a, b) ; print_matrix(c, 2, 2) ; return 0 ;
}
void add_matrix(double c[][MAX], double a[][MAX], double b[][MAX]) {
int i, j ; for(i=0;i<MAX;i++)
for(j=0;j<MAX;j++)
c[i][j] = a[i][j] + b[i][j] ; return ;
}
void print_matrix(double a[][MAX], int rows, int cols) {
int i, j;
for(i=0;i<rows;i++) { fprintf(stdout,"[") ;
for(j=0;j<cols;j++) fprintf(stdout,"\t%f", a[i][j]) ; fprintf(stdout,"]\n") ;
}
☆ double 型の配列で表現する
/* $Id: matrix_1.c,v 1.2 2006-06-15 19:57:36+09 naito Exp $ */
#include <stdio.h>
#define MAX (10)
void print_matrix(double *, int, int) ;
void add_matrix(double *, double *, double *, int, int) ; int main(int argc, char **argv)
{
double a[MAX*MAX] = {1.0, 2.0, 2.0, 3.0} ; double b[MAX*MAX] = {1.0, 2.0,
2.0, 3.0} ; double c[MAX*MAX] ;
print_matrix(a, 2, 2) ; print_matrix(b, 2, 2) ; add_matrix(c, a, b, 2, 2) ; print_matrix(c, 2, 2) ; return 0 ;
}
void add_matrix(double *c, double *a, double *b, int rows, int cols) {
int i, j ;
for(i=0;i<rows;i++) for(j=0;j<cols;j++)
c[i*cols+j] = a[i*cols+j] + b[i*cols+j] ; return ;
}
void print_matrix(double *a, int rows, int cols) {
int i, j;
for(i=0;i<rows;i++) { fprintf(stdout,"[") ;
for(j=0;j<cols;j++) fprintf(stdout, "\t%f", a[i*cols+j]) ; fprintf(stdout,"]\n") ;
}
fprintf(stdout,"\n") ; return ;
}
☆ 構造体を使ったもの
/* $Id: matrix_2.c,v 1.2 2006-06-15 20:00:08+09 naito Exp $ */
#include <stdio.h>
#define MAX (10) struct matrix {
int rows ; int cols ;
double data[MAX*MAX] ; } ;
typedef struct matrix matrix ; void print_matrix(matrix) ;
void add_matrix(matrix *, matrix, matrix) ; int main(int argc, char **argv)
{
matrix a = {2, 2, {1.0, 2.0,
2.0, 3.0}} ; matrix b = {2, 2,
{1.0, 2.0, 2.0, 3.0}} ;
matrix e ;
print_matrix(a) ; print_matrix(b) ; add_matrix(&c, a, b) ; print_matrix(c) ; return 0 ;
}
void add_matrix(matrix *c, matrix a, matrix b) {
int i, j ;
c->rows = a.rows ; c->cols = a.cols ;
for(i=0;i<c->rows;i++) for(j=0;j<c->cols;j++)
c->data[i*c->cols+j] = a.data[i*a.cols+j] + b.data[i*b.cols+j] ; return ;
}
void print_matrix(matrix a) {
int i, j;
for(i=0;i<a.rows;i++) { fprintf(stdout,"[") ;
for(j=0;j<a.cols;j++) fprintf(stdout, "\t%f", a.data[i*a.cols+j]) ; fprintf(stdout,"]\n") ;
}
fprintf(stdout,"\n") ;
★ 基本変形行列
• Di(λ): aii=λ,ajj = 1 (j6=i) Di(λ)−1=Di(λ−1), detDi(λ) =λ.
i 行目(i列目)をλ倍する.
• Tij: akk= 1 (k6=i, j),aii =ajj = 0, aij=aji= 1.
Tij−1=Tij, detTij= (−1)|i−j|.
i 行目とj 行目(i行目とj 行目)を入れ替える.
• Eij(λ): akk= 1,aij =λ.
Eij(λ)−1 =Eij(−λ), detEij(λ) = 1.
j 行目のλ倍をi行目(i列目のλ倍をj 列目)に加える.
Di(λ) =
i
<
1 0 . . . 0
0 . .. ...
... . .. ...
0 . . . λ . . . 0
... . .. ...
... . .. 0
0 . . . 0 1
< i
Tij =
i
<
j
<
1 . . . 0 ... . .. . . ... 0 · · · 0 · · · 1 · · · 0
... . .. ...
0 · · · 1 · · · 0 · · · 0 ... . . .. ...
0 . . . 1
< i
< j
Eij(λ) =
i
<
j
<
1 . . . 0 ... . .. . . ... 0 · · · 1 · · · λ · · · 0
... . .. ...
0 · · · 0 · · · 1 · · · 0 ... . . .. ...
0 . . . 1
< i
< j
● 実習内容
1. 行列のスカラー倍・積を計算するプログラムを書きなさい.
2. Gauss-Jordanの消去法,およびGauss の消去法を用いて,以下の連立方程式を解くプログラムを書
きなさい.
3x+12y+9z=3, 2x+ 5y+4z=4,
−x+ 3y+2z=5.
2x+2y+ z=3, x+ y+2z=4, y+ z=5.
3.
A=
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
, b=
sin(π/(n+ 1)) sin(2π/(n+ 1))
...
sin((n−1)π/(n+ 1)) sin(nπ/(n+ 1))
に対して,連立一次方程式Ax=bの解を求めるプログラムを書きなさい. （たとえば,n= 10とか でよい.）
4. Gauss-Jordanの消去法, Gauss の消去法で, 行列が正方行列でない場合や full rankでない場合に, 可解性を判定したり,可解であるときには,解を一組求めたり,核の基底をもとめることができるよう にアルゴリズムを改良しなさい.
5. Gauss-Jordanの消去法を用いて,以下の連立方程式が可解かどうかを判定し,可解であれば,その解
をひとつ求めなさい. また, その連立方程式をあらわす行列が生息でない場合, その核の基底を一組 求めなさい.
2 2 1 1 1 2 0 0 1
x y z
=
−4 1 2
2 2 1 1 1 2 0 0 1
x y z
=
2 3 4
1 1 1 1 1 1 1 1 1
x y z
=
1 1 1
0 0 0 0 0 0 0 0 0
x y z
=
2 3 4
1 2 1
1 1 4
!
x y z
= 3 4
!
1 2
−3 2
−1 6
x y
!
=
3 4 5
1 2
−3 2
−2 4
x y
!
=
1 1 2
