計算数学基礎（2011 後期）

シラバス 00-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

この講義について作成日 : October 6, 2011 Version : 1.1

担当教員：
川平 友規（かわひら ともき，講義担当，A441 号室，kawahira@math.nagoya-u.ac.jp）
佐藤 猛（さとう たけし，演習担当，A359 号室，sato@math.nagoya-u.ac.jp）
講義のウェブサイト：
http://www.math.nagoya-u.ac.jp/~kawahira/courses/11W-keisan.html
配布されたプリントが pdf 形式でダウンロードできます．また，毎週の進捗状況についてコメン
トしていきます．
本講義の目的（コースデザインより）： 本講義の目的は，数理科学の問題に対してコンピュータを活用するための基礎知識を習得することである．具体的には，数式処理ソフトウェア Mathematica
を用いて，数理科学の諸問題に取り組む．
講義日と授業内容（予定）
：
第1回
第2回
第3回
第4回
第5回

10/6
10/13
10/20
10/27
11/10

基礎編（Mathematica 入門）：
端末の設定，代数（記号）演算，方程式，
数列，リストの操作，微分・積分，
線形代数，平面グラフィックス，空間グラフィックス

第6回
第7回
第8回
第9回
第 10 回

11/17
11/24
12/1
12/8
12/15

応用編（「可視化」をテーマに）：テイラー展開，フーリエ展開，
ベクトル場，等位面，微分方程式，
複素関数，乱数，ブラウン運動，
パーコレーション，素数とゼータ

第 11 回
第 12 回
第 13 回
第 14 回

12/22
1/12
1/19
1/26

発展編（プログラミングの初歩)：
計算誤差の問題，ループ（Do, For, While）
力学系，カオス，フラクタル
連分数，セル・オートマトン

講義の進め方： 1 限は教室での講義，2 限は下記の端末室での実習となります．

• 講義：理 1 号館 109 号室• 演習：工学部 7 号館 4 階情報メディア教育センター（主センター）
・端末室 B
1 限の最初に出席確認を兼ねたクイズ（小テスト）を行います．また 2 限の実習中に，その日の課題をメールに添付し提出（送信）してください．
成績の評価方法： 1・2 限の出席とクイズ，課題を総合的に判断して評価します．
（1 限の遅刻は 3
回目以降欠席とカウント．欠席 5 回で単位放棄とみなします．
）
課題の提出方法： 演習で出される課題は，次の要領でメールの添付ファイルとして送信してくだ
さい．
● 送信先：comp0-2011@math.nagoya-u.ac.jp （comp0 の 0 はゼロ．川平と佐藤に届きます．
）● 締切：2 限の講義終了時．その後は一切受け付けません．
● メールの様式：メールの件名と添付ファイル名はともに半角英数字で「講義の回数-学籍番号」
URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

シラバス 00-2

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jpの形とします．たとえば，第 5 回目の課題を学籍番号 060123456 の人が提出するときは
件名： 05-060123456
添付ファイル名： 05-060123456.m （ファイル形式については後述）のようにすること．
！
！注意！
！ 不正行為を防ぐために，学生間のコンピュータファイルの受け渡しは禁止します．
その他，不正行為には厳しく対処するので注意すること．
その他：(1) 学内のコンピュータの利用法に関する情報は情報メディア教育システムのページ
http://www.media.nagoya-u.ac.jp/ にまとめられています．このページは「ブックマーク」に登録しておくとよいでしょう．
(2) Mathematica を自習・活用したいとき：平日の 8:45—17:00 の間，各端末室のパソコンは講義に使われている時間を除いて自由に利用することができます．たとえば今学期，上記の時間帯であれば主センターの端末室 A・B のいずれかが利用できる予定です．詳しくはhttp://www.media.nagoya-u.ac.jp/lectures/ の講義予定表を参照してください．他の端末室の利用スケジュールも確認できます．
参考書：• 日本 Mathematica ユーザー会, 「入門 Mathematica」（東京電機大学出版局）• 榊原進, 「はやわかり Mathematica」（共立出版）
オフィスアワー： 授業中・授業前後の質問は大歓迎です．それ以外の時間に質問したい場合は，ぜひオフィスアワー（教員ごとの質問受付時間）を活用してください．私のオフィスアワーは，Cafe
David（カフェ・ダヴィド）という合同のオフィスアワー内に設定しています．Cafe David は月曜から金曜の 12:00-13:30 にオープンし，コーヒー・紅茶を無料で提供しています．私の担当は月曜
日です．その他のスタッフ（名札をつけています）も待機しているので，自由に質問していただ
いてかまいません．数理学科の学部生・院生もたくさんいるので，勉強のこと，進路のことなど，情報収集の場としても活用してください．
謝辞： 授業を準備するにあたって，2008-2009 年度にこの講義を担当された糸健太郎さん，2010年度に担当された粟田英資さんにはたいへん貴重なアドバイスをいただきました．とくにこのシ
ラバスは糸さんが授業で使われたものをベースにしております．また内藤久資さんにはシステム関連でご協力いただきました．みなさん，どうもありがとうございました．
よく使う記号など：数の集合(1) C: 複素数全体(4) Z: 整数全体

(2) R: 実数全体
(5) N: 自然数全体(3) Q: 有理数全体(6) x ∈ R: x は実数

ギリシャ文字(1) α: アルファ (2) β: ベータ (3) γ, Γ: ガンマ (4) δ, ∆: デルタ (5) ϵ: イプシロン(6) ζ: ゼータ(7) η: エータ (8) θ, Θ: シータ (9) ι: イオタ(10) κ: カッパ(11) λ, Λ: ラムダ (12) µ: ミュー (13) ν: ニュー(14) ξ, Ξ: クシー (15) o: オミクロン(16) π, Π: パイ (17) ρ: ロー(18) σ, Σ: シグマ (19) τ : タウ(20) υ, Υ: ウプシロン(21) ϕ, Φ: ファイ (22) χ: カイ(23) ψ, Ψ: プサイ (24) ω, Ω: オメガ

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 01-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

Mathematica 入門
作成日 : October 6, 2011 Version : 1.3

今日やること(1) Mac にログインする．
Mac が英語環境になっているときは，次のようにして日本語環境に変更する：
メニューバーの「アップルマーク」→[System Preference...]→[International] この中の Languages の中の [日本語] を [English] の上に移動させる．このウィンドウを閉じて，ログアウ
トし，再度ログインする．

(2) メールソフト Thunderbird の設定をする．
(a) まず Dock 1 のアイコンをクリックして Thunderbird を起動する．アイコンが Dock に無い場合は，「TEMPHDD」→「アプリケーション」→Thunderbird.app と辿る．このアイコンは Dock の中にドラッグ＆ドロップしておくとよい(b) Thunderbird でメールの送受信設定をする．設定の方法は別紙参照のこと．
(c) 設定が終わったら，まず自分にテストメールを送ってみる．
(d) 自分へのメールが無事に送信・受信できること確認したら，次にcomp0-2011@math.nagoya-u.ac.jpにテストメールを送る．メールの件名は「00-060123456」（半角で．00-のあとは自分の学籍番号にかえる）の形とし，本文に名前（ふりがな）を書く．
(3) Mathematica の起動，設定，実行．
(a) Mathematica を起動する．Mathematica のアプリケーション本体は Thunderbird と同じところにある．
（このアイコンも Dock の中にドラッグ＆ドロップしておくとよい．
）(b) Mathematica が日本語環境でない場合：メニューバーの「Mathematica」→「Preferences」→ 「Interface」で，「Language for menus and dialog boxes」にあるドロップメ
ニューから「Japanese」を選ぶ．さらに Mathematica を一度終了し，再度起動する．
(c) Mathematica の文字（フォント）を適宜大きくする（これは暇なときにやればよい）：
メニューバーの「Mathematica」→「環境設定．．
」→ 「詳細」→ 「オプションイン
スペクタを開く」→ 「書式設定」→ 「フォント設定」→ 「フォントサイズ」これを20 ∼ 30 ぐらいの数にする．
(d) 次ページの課題をがんばる．
(e) Mathematica の「名称未定義-1」というファイルをメニューバーの「ファイル」→「別名で保存...」として，01-060123456.nb のような名前で保存する．保存場所は「TEMPHDD」→ （家のアイコン）→「書類」の中とする．必要に応じて新規フォルダを作ってもよい．
（Mathematica はまだ閉じちゃだめ．
）(f) 次に同じファイルを今度はパッケージ形式で保存する：すなわち，「別名で保存...」から「フォーマット」のタブで Mathematica パッケージ (*.m) を選び，01-060123456.m という名前にして書類フォルダに保存する．

1

画面下のアイコンが並んでいるところ

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 01-2

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

課題を提出する際の注意！
！通常はノートブック形式（*.nb）で保存するが，このままだと
ファイルサイズが非常に大きいので，課題としてメールに添付する際にはパッケージ形式（*.m）で保存してこちらを送ること．

(g) 課題の提出：この 01-060123456.m をメールに添付し comp0-2011@math.nagoya-u.ac.jpに送る．ただしメールの件名は「01-060123456」とすること．
今日の課題 (10/6).
右側の式を 1 行ずつ入力・評価（evaluate）しノートブックを作成せよ．
式の評価は Shift + Return （ Shift を押しながら Return ）である．
和：
積と商：
分数：
べき乗：
四則と定数
階乗：
定数 π ：
π の数値化：π を 100 桁：
1000 番目の素数：

1 + 1
15*(10 + 3)/12
1/3 + 1/12
3^10
100!
Pi
N[Pi]
N[Pi, 100]
Prime[1000](a + b)10 の展開，結果を q とおく ：q をもう一度因数分解 ：方程式を解く：
文字式と方程式
文字が入っても OK ：1 の 10 乗根：同じ方程式を数値的に解く：
sin x を x で微分：
cos xy 2 を y で偏微分：
微分・積分
cos xy 2 を x で偏微分：∫1
定積分 0 x2 dx ：q = Expand[(a + b)^10]
Factor[q]Solve[x^2 - x - 2 == 0, x]
Solve[x^2 - x + c == 0, x]
Solve[x^10 == 1, x]
NSolve[x^10 == 1, x]

D[Sin[x], x]D[Cos[x*y^2], y]D[Cos[x*y^2], x]Integrate[x^2, {x, 0, 1}]

sin x の −π ≤ x ≤ π でのグラフ： Plot[Sin[x], {x, -Pi, Pi}]sin(ex ) の −5 ≤ x ≤ 5 でのグラフ: Plot[ Sin[Exp[x]], {x, -5, 5}]√グラフ
sin x2 + y 2 の −10 ≤ x, y ≤ 10 でのグラフ：
Plot3D[ Sin[Sqrt[x^2 + y^2]], {x, -10, 10}, {y, -10, 10}]
研究• 課題が提出できたら，ノートブックをさらに別名で保存し，数値や関数を変えてどんな結果がでるか試してみよ．
• ヘルプ −→ ドキュメントセンターと辿ってドキュメントセンターを開き，Mathematica の機能についていろいろ勉強せよ．
• ドキュメントセンターの検索ウインドウから N Plot Sin など，今日使ったコマンドを検索してみよ．
URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 02-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

式の変形と値
作成日 : October 13, 2011 Version : 1.2

役立つキー操作
やり直し：
コピー：
貼り付け：
切り取り：

Command + Z
Command + C
Command + V
Command + X

式の評価（実行）：
前の入力をコピー ：単語を補完：√
根号の入力：

Shift + Return
Command + L
Command + K
Command + 2

注．Windows や UNIX の場合 Command は Ctrl となる．

入力時の注意(1) 必ず半角英数字で入力．
（全角漢字も使えるが混乱のもとになる．
）](2) 入出力を消したいときはノートブック右側のをマウス等で選択し， Delete .
(3) 数の掛け算は x*y でも x y（間に半角空白をはさむ）でもよい．
(4) 計算の優先順位を表す括弧は常に ( )：1/{1 + (2 + 3)} ではなく，1/(1 + (2 + 3)). {}はリストやベクトル（{a,b,c}など）に，[] は関数（Sin[x] など）に用いられる．
(5) Mathematica 内ですでに定義されている関数（Sin[x], Exp[x] など）は大文字で始まる．
(6) Mathematica は一度代入された値をずっと覚えている．もし a = 1; b = a + 1 と入力・評価されたら，そのあと a と b の値はそれぞれ 1 と 2 であり続ける．再度別の意味で使うときは，かならず値をクリアしよう：Clear[a, b] もしくは a = .; b = .
(7) 直前の出力は % で参照できるが，混乱のもとになるので課題には使わないこと．
基本的な関数と定数（関数）注．下記の定数は「厳密値」として記号のまま扱われる．近似値がほしいときは N 関数などで適宜数値化する．

円周率（定数）：
自然対数の底（定数））：
虚数単位（定数）：
黄金比（定数）
：
平方根：
三角関数：

逆三角関数：

Pi
E
I
GoldenRatio
Sqrt[ ]
Sin[ ]
Cos[ ]Tan[ ]
ArcSin[ ] など

指数関数：
自然対数 log a：
対数 logb a ：数値化：
m 桁で数値化 ：x 以下の最大の整数 ：a, b, c, . . . の最大値 ：a, b, c, . . . の最小値 ：Exp[ ]
Log[a]
Log[b , a]N[ ]N[ , m]Floor[ ]Max[a,b,c,...]
Min[a,b,c,...]

注意. 以下とくに断らないかぎり，「近似値」といったら N 関数による値とする．
問題 1. Mathematica を用いて計算せよ．

(1) π e , eπ , π + 20 の近似値をそれぞれ有効数字 7 桁で求め，もっとも小さなものを求めよ．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 02-2

担当教員：川平 友規 研究室：A441√

Email: kawahira@math.nagoya-u.ac.jp√√√
1√ , d = 3 とする．これらの近似値をもとめ，3, b = 15 − 2, ．c = √5+ 3
最大のものと最小のものを決定せよ．
（神戸薬大・改）(2) a =√√√a b
2i, b = 3 − 2i のとき，c = + と d = a4 − 2a2 それぞれの厳密値および
b a
有効数字 20 桁の近似値を求めよ．
（’98 センター追試・改）√√p(4) 42 + 12 6 の整数部分を p，小数部分を q とするとき，p, q, r =の値を厳密値q(q + 4)および近似値で求めよ．
（成蹊大・改）(3) a =√

2+

3+

文字式の操作問題 2. 順に入力して結果を確かめよ．
整数係数での因数分解：式の展開，結果を s とせよ：
s を x でまとめる：ｓ に x = 2 を代入：
ｓ に x = 2, y = −3 を代入：
1t0 = 3
とおく：
x −1
t0 を部分分数に分解，結果を t1 とおく：
t1 を通分，結果を t2 とおく：
t2 の式の分母だけ展開する：
t2 の式を簡単にする：
Simplify で簡単に（単なる計算）：
Simplify で簡単に（公式を適用）：
FullSimplify はもっとすごい：

Factor[x^10-1]
s = Expand[(1 + x + y)^3]
Collect[s, x]
s/.{x -> 2}s/.{x -> 2, y -> -3}
t0 = 1/(x^3-1)
t1 = Apart[t0]t2 = Together[t1]ExpandDenominator[t2]Simplify[t2]Simplify[Sqrt[512]/4]
Simplify[Sin[x]^2 + Cos[x]^2]
FullSimplify[Cos[x] + I Sin[x]]問題 3. Mathematica を用いて計算せよ．

(1) 問題 1 の (3)(4) の入出力部分にもどり，c, d, q, r の厳密値に Simplify（もしくは FullSimplify）
を施せ．
(2) x4 − 3x2 y 2 + y 4 を因数分解せよ．その結果に，x = −1 を代入せよ． （04 名古屋経大・改）√√(3) (x + y)(y + z)(z + x) + xyz を展開し，さらに因数分解せよ．その式に x = 2 および y = 3
を代入せよ．
（00 名城大・改）
(4)

1
a
b
c
=
+
+
が恒等式であるとき，a, b, c の値を求め(x + 1)(x + 2)(x + 3)x+1 x+2 x+3よ．
（97 静岡理工科大）

(5)

1
2x
3x2
a
b
c
+
+
=
+
+
が恒等式であるとき，a, b, c
x + 2 (x + 2)2 (x + 2)3
x + 2 (x + 2)2 (x + 2)3の値を求めよ．
（92 久留米工大）

(6)

b2
c2
a2
+
+
を簡単にせよ．
(a − b)(a − c) (b − c)(b − a) (c − a)(c − b)

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/（02 横浜市大・改）名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）
担当教員：川平 友規 研究室：A441

演習問題 02-3
Email: kawahira@math.nagoya-u.ac.jp

課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイル名を 02-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：02-060123456
本文：もし今日の感想・意見・要望などがあれば書いてください．

研究(1) 複素関数論で習った下記の複素数の値と右での Mathematica による出力結果（とその近似値）を比べよ：√−1： Sqrt[-1](−1)1/3 ： (-1)^(1/3)
eπi/4 ： Exp[Pi*I/4]
sin i ： Sin[I]arcsin 3 ： ArcSin[3]log i ： Log[I]log(1 + i) ： Log[1 + I]
(1 − i)i ： (1-I)^I
√ √
√
6+1
6+3
√
(2) FullSimplify を用いて a =+√を簡約化しても，2 重根号が残ってしま
6−1
6+1
√
うことを確認せよ．次に 5a を簡約化するとどうなるか？
（97 名古屋学院大・改）(3) 次の問題をMathematica
が厳密値で解けるか試してみよ．近似値の場合はどうか？：x =
√
√
√
√
3
3
7 + 5 2，y =7 − 5 2 のとき，x + y =, xy =となるから，x4 + y 4 =, x5 + y 5 =である．
√(Hint: y 3 = 7 − 5 2 < 0 であるから，y < 0 としたいのだが…？)（84 早大・改）(4) ドキュメントセンターをひらき，Expand や Together に関連してどのような組み込み関数があるか調べよ．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 03-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

リストと数列作成日 : October 20, 2011 Version : 1.2

Tips : 評価の放棄
Mathematica に計算させていると，予想外に時間がかかってしまうことがある．入力式の評価を放棄 (abort) させるには， Command + . （ピリオド）を用いる．
（Windows では Alt + . となる． メニューから「評価」→ 「評価を放棄」を選んでもよい．
）
リストとは？
リストとはベクトル，行列，テンソル，数列，集合などなど，さまざまなものを表現できる
Mathematica 独自のデータ構造（データのいれもの）であり，中括弧 { } を用いて

{a, b, c}{{a, b}, {c, d}},{Sin[x], Cos[x], Tan[x]}のように表される．
リストの計算リストにはさまざまな演算が大胆に適用できる．
リスト v1 を定義：
リスト v2 を定義：
ベクトル的な和：
ベクトル的な定数倍：
後ろからの定数倍：
各要素に 1 を足す:
文字 x から各要素を引く:

v1 = {a, b, c}
v2 = {p, q, r}
v1 + v2
100*v1v1*100v1 + 1x - v1

各要素を 3 乗：
5 を「リスト乗」：
指数関数を各要素に作用：リストの「リスト乗」：
要素ごとの積：
要素ごとの商：ベクトルとしての内積：

v1^3
5^v1
Exp[v1]
v1^v2
v1*v2
v1/v2
v1.v2

問題 1. 上の表のコマンドを入力し，どのような結果が得られるか確認せよ．
問題 2. u = {1, 2, 3, 4, 5}と定義する．これと上の表のコマンドを活用して，{}{}
x2 x3 x4 x52
345x − 1, (x − 1) , (x − 1) , (x − 1) , (x − 1) ,x,,,,と等しいリストをそ2! 3! 4! 5!
れぞれ作成せよ．ただし，n の階乗は n!もしくは Factorial[n] で計算できる．
リストの生成リスト {1, 2, 3, 4, 5} を定義：
リスト {4, 5, . . . , 10} を定義：
4 ≤ x ≤ 10 で 0.7 刻みのリスト：
一般項でリストを作る（２乗）
：
（文字式）
：
リストの結合：リストの和集合：
リストの共通部分：

v1 = Range[5]v2 = Range[4,10]v3 = Range[4,10, 0.7]
v4 = Table[n^2, {n,1,10}]v5 = Table[x^n, {n,0,10}]
Join[v1, v2]
Union[v1, v2]
Intersection[v1, v2]

問題 3. 上の表のコマンドを入力し，どのような結果が得られるか確認せよ．
()
1 n
問題 4. Table を用いて，自然対数の底 e への収束列 an = 1 +の第 10 項目まで近似値（Nn関数を施した値）からなるリスト an を作れ．
URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 03-2

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

問題 5. Table を用いて，x − 1, x2 − 1, · · · , x10 − 1 からなるリストと，その各要素を因数分解した結果からなるリストを作れ．
リストのリストのリストの…
リストのリストを定義：
行列 (matrix) として表示：
表 (table) として表示：
一般項でリストを作る（和）
：
（文字式）
：

m1 = {{a, b}, {c, d}, {f, g}}
m1 //MatrixForm
もしくは MatrixForm[m1]m1 //TableForm
もしくは TableForm[m1]m2 = Table[i + j, {i,1,5}, {j,1,5}]m3 = Table[x^i*y^j, {i,1,5}, {j,1,5}]

問題 6. 上の表のコマンドを入力し，どのような結果が得られるか確認せよ．
問題 7. (九九の表) リストの TableForm を用いて，1 から 9 までの積の表を作れ．また，1 から21 までの奇数同士の積を表にせよ．
問題 8. (一覧表の作成) 「{n, n^2, n^3} の形のリスト」からなるリスト
pow = Table[{n, n^2, n^3}, {n,1,10}]を作成せよ．さらに，これを MatrixForm および TableForm を用いて表せ．
リストから要素を取り出す
再度 2 乗 (square) の列：
リスト sq の 7 番目の要素：
再度，和 (wa) の行列：
再度行列として表示：
行列 wa の 2 行 3 列目：

sq = Table[n^2, {n,1,10}]
sq[[7]]
wa = Table[i + j, {i,1,5}, {j,1,5}]
wa //MatrixForm
wa[[2,3]]

問題 9. 上の表のコマンドを入力し，どのような結果が得られるか確認せよ．
問題 10. (収束列と誤差の一覧表)(1) 問題 4 で作ったリスト an を用いて，an と e との誤差 |an − e| からなるリスト anError を作れ．ただし，数 x の絶対値は関数 Abs[x] で与えられる（Hint: たとえば |a5 − e| はリ
スト an を使うと Abs[an[[5]]- E] で表される．
）(2) 問題 8 を応用して, an と誤差 |an − e| を比較する一覧表を作成せよ．
数列の有限和と無限和
またまた 2 乗の列：
1 から 10 までの 2 乗の和：
一般項：
1 から 10 までの −2 乗の和：
−2 乗の無限和（級数）：
10 次多項式を作る：
じゃあこの級数は？
：数列の積は Product を使う：
またまた和の行列：：
これらすべての和：
一般項：

sq = Table[n^2, {n, 1, 10}]s1 = Sum[n^2, {n, 1, 10}]s2 = Sum[k^2, {k, 1, n}]s3 = Sum[1/n^2, {n, 1, 10}]s4 = Sum[1/n^2, {n, 1, Infinity}]s5 = Sum[x^n/n!, {n, 0, 10}]s6 = Sum[x^n/n!, {n, 0, Infinity}]p = Product[2*n - 1, {k, 1, 10}]wa = Table[i + j, {i,1,5}, {j,1,5}]
wa //TableForm
Sum[i + j, {i,1,5}, {j,1,5}]
Sum[i + j, {i,1,n}, {j,1,m}]

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 03-3

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

問題 11. 上の表のコマンドを入力し，どのような結果が得られるか確認せよ．

1 1
1+ +· · ·+1! 2!
n!
の第 10 項目までの近似値からなるリストを作成せよ．また，問題 10 と同様にして，an , |an −e|, bn , |bn − e| を比較する一覧表を作成せよ．
問題 12. (収束速度の比較) Table と Sum を組み合わせて，e への収束列 bn = 1+

問題 13. (大学入試問題) Mathematica を用いて計算せよ．

(1) 次の和を求めよ：(2) 数列100∑

3√. （近似値だけでよい）√
k+ k−1k=2（99 兵庫大）22
2,,, . . . の第 n 項のまでの和を求めよ．
22 − 1 32 − 1 42 − 1（00 日本福祉大）(3) 連続する m 個の奇数 1, 3, 5, . . . , 2m − 1 の中から，異なる 2 つの数をとって積をつくる．
こうして得られる m C2 通りの積すべての和を求めよ．
（99 浜松医大）(4) 初項 4，公比 8 の等比数列の第 n 項までの和を Sn で表す．このとき，S99 は カキ けたの
整数，S100 は クケ けたの整数である．さらに，S1 から S100 の間で， クケ 以下のすべてのけた数の Sn が存在するか確認せよ．また，Sn が初めて 150 けたを超えるのはいつかを，試行錯誤で確定せよ．
（92 センター本試・改）(Hint: Sn の桁数は 1 + log10 Sn の整数部分で与えられる．その一覧表を作れ．)
課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイル名を 03-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：03-060123456
本文：もし今日の感想・意見・要望などがあれば書いてください．
おねがい：退出時，iMac をログアウトでなくシャットダウンしてください．

研究(1) i, j, k が 1 から
∑n までの値をすべて動くとき，i + j + k の和を求めよ．
すなわち，(i + j + k) を計算せよ．
1≤i,j,k≤n(2) 一般項 pn = n2 + n + 41 で与えられる数列が初めて素数でなくなる n を求めよ．
（ある mが素数かどうかを判定する組み込み関数は PrimeQ[m] である．
）
(3)
sn :=

n
∑
k=1

cos

2πk
,
n

pn :=

n
∏
k=1

cos

2πk
n

とおく．n = 5, 6, 7, 8, 9 について sn と pn を計算せよ．また FullSimplify によって，snと pn の一般項を予測せよ．
（じつは，Mathematica は一般項を知っている．
）(4) ドキュメントセンターをひらき，リストに関連してどのような組み込み関数があるか調べよ．
URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 04-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

方程式とグラフ
作成日 : October 27, 2011 Version : 1.1

Tips : 入力結果を省略させるには？長い式やリスト（数列や行列など）を Table などで定義すると，入力結果の表示だけで画面がいっぱいになってしまう．このような表示を省略したければ，入力式の最後にセミコロン ; をつけるとよい．

2 次方程式とグラフ
方程式は組み込み関数 Solve で解ける．答えは代入形（置換則リスト）で返ってくるので，扱いには注意が必要である．また解の様子をグラフで吟味するには，Plot を使う．

x2 − x − 1 = 0 を解き，結果を sol とおく：
厳密解から近似値へ：
いきなり数値解（近似解）を求める：
有効数字を 100 桁に指定：グラフで確認（範囲は −1 ≤ x ≤ 2）
：

[1]
[2]
[3]
[4]
[5]

sol = Solve[x^2 - x - 1 == 0, x]N[sol]
NSolve[x^2 - x - 1 == 0, x]
NSolve[x^2 - x - 1 == 0, x, 100]
Plot[x^2 - x - 1, {x, -1, 2}]

問題 1. 上の表のコマンドを入力し，結果を確認せよ．
問題 2. (解の図示) Solve 関数および NSolve 関数を用いて方程式 x3 − 19x + 30 = 0 を解け．また f (x) = x3 − 19x + 30 のグラフを描け．ただし，x 軸との交点が分かるように適当な描画範囲を選ぶこと．
方程式の解を活用する
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]

代入（置換則）の復習．x3 に x = 5 を代入 ：さっきの方程式 x2 − x − 1 = 0 の解：解だけからなるリストを作成：
解のリストに直接代入：
値を確認：
：
3
3
式 t = x − 1/x の値を求めよう：
t に解を両方とも（a と b）代入：
値を整理：
t に sol のひとつ目の解（a だけ）を代入：
ふたつ目の解（b だけ）を代入：

x^3 /. {x -> 5}
sol
x /. sol{a, b} = x /. solabt = x^3 - 1/x^3
ans = t /. sol
Simplify[ans]ans1 = t /. sol[[1]]ans2 = t /. sol[[2]]問題 3. 上の表のコマンドを入力し，結果を確認せよ．
問題 4. 方程式 x2 − x − 5 = 0 の解を α, β (α < β) とするとき，α4 + β 4 および α4 − α2 − 30 を計算せよ．([9]-[11] を使う．)
問題 5. (大学入試問題) 結果には必要に応じて Simplify もしくは FullSimplify を施すこと．
11(1) x2 − 3x + 1 = 0 のとき，x3 + 3 , x5 + 5 の値を求めよ．
（97 帝塚山大）x
x
1(2) x2 + 2 = 4 のとき，次の値を求めよ．ただし，1 < x < 2 とする．
x111（ア）x2 − 2（イ）x3 − 3（ウ）x5 − 5
（01 大同工大）
xxx(Hint: sol = Solve[...] で解いたあと，sol に N 関数を施す．1 < x < 2 となる解が i 番目であれば，（ア）〜（ウ）に sol[[i]] を代入．)

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 04-2

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

連立方程式
[17]
[18]
[19]

[20]
[21]

直線 vs 直線：
y = の形のグラフで見る：
陰関数でグラフを描く：（範囲は −5 ≤ x, y ≤ 5）：
直線 vs 放物線：
陰関数のグラフ：（範囲は −4 ≤ x, y ≤ 4）：

s1 = Solve[{x + y == 1, x - y == -2}, {x, y}]
g1 = Plot[{-x + 1, x + 2}, {x, -5, 5}]
g12 = ContourPlot[{x + y == 1, x - y == -2},
{x, -5, 5}, {y, -5, 5}]s2 = Solve[{x + y == 1, y == x^2 - 1}, {x, y}]g2 = ContourPlot[{x + y == 1, y == x^2 - 1},{x, -4, 4}, {y, -4, 4}]

問題 6. 上の表のコマンドを入力し，結果を確認せよ．
問題 7. (直線 vs 円) [20][21] を応用して，連立方程式 x + y = 1, x2 + y 2 = 2 を解き，陰関数としてのグラフを −3 ≤ x, y ≤ 3 の範囲で描け．
√
3
1
問題 8. (連立三角方程式) 同様にして連立方程式 sin x + sin y =, cos x + cos y = を解き．
22陰関数としてのグラフを −2π ≤ x, y ≤ 2π の範囲で描け．
（63 滋賀大・改）連立方程式の解の活用
[22]
[23]
[24]
[25]
[26]
[27]
[28]

さっきの連立方程式（直線 vs 放物線）の解：
解をベクトル v1 と v2 に代入：
値を確認：：
解を用いた計算：
両方の解（v1 と v2）を代入：
最初の解（v1）だけ代入：s = s2{v1, v2} = {x, y} /. sv1v2
h = x^5 + y^5;ans = h /. s
ans = h /. s[[1]]問題 9. 上の表のコマンドを入力し，どのような結果が得られるか確認せよ．
問題 10. (3 元連立方程式)（明治薬大・改）(1) x, y, z が方程式 x+y +z = 2, x2 +y 2 +z 2 = 6, x3 +y 3 +z 3 = 8 を満たすとき，xy +yz +zx
および xyz の値をもとめよ．さらに，x2 − y + z が最大になるとき，x, y, z の値を求めよ．
(2) 下のコマンドを入力して 3 つの方程式の陰関数グラフ（曲面）を描き，これらの曲面の共通交点が６つあることを確かめよ．：
gr = ContourPlot3D[{x + y + z == 2, x^2 + y^2 + z^2 == 6, x^3 + y^3 + z^3 == 8},{x, -3, 3}, {y, -3, 3}, {z, -3, 3},ContourStyle -> {Yellow, Green, Orange},
Mesh -> None, AxesLabel -> Automatic ]課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイ
ル名を 04-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：04-060123456
本文：もし今日の感想・意見・要望などがあれば書いてください．
おねがい：退出時，iMac をログアウトでなくシャットダウンしてください．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 05-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

微分・積分・いろんなグラフ
作成日 : November 10, 2011 Version : 1.2]
Tips : ノートブックはセル（）で分割・グループ化されている．タイトル・小見出しなどのセルを作るには，以下のショートカットがおすすめ：
でかい字のタイトル (title)：
節 (section），横線が入る：
小見出し ((sub)subsection)：

Command + 1
Command + 4
Command + 5 もしくは 6◎今日の課題◎：問題の前にある [1][2]…等のコマンドをまず入力し，結果を確かめた上で問題にとりかかること．

関数の定義（即時割り当て，自分で関数を作る）
[1]
[2]
[3]
[4]
[5]

関数の定義（即時割り当て）：
複素数 x = 1 + i を代入：
リスト {2,3} を代入：
2 変数関数：値を代入：

f[x_] = x^2 + 1
f[1 + I]f[{2, 3}]g[x_, y_] = x y + 1
g[3, 3]

1 変数関数の微分
[6]
[7]
[8]
[9]
[10]
[11]
[12]

関数 f を再定義：
x で微分（文字式 d1 とおく）：
x = 2 での微分の値：
2 階微分：
こららの関数のグラフを描く：
f の導関数を定義：
x = 2 での導関数の値：

f[x_] = x^3 - 3 x;
d1 = D[f[x], x]
d1 /. {x -> 2}
d2 = D[f[x], {x, 2}]
Plot[{f[x], d1, d2}, {x, -4, 4}]g[x_] = D[f[x], x]g[2]問題 1. (極値を求める) 上の [11] で定義した f[x] の導関数 g[x] が 0 となる x を Solve 関数
を用いて求めよ．さらに，そのときの関数 f[x] の値を求めよ．(Hint: 前回のプリント [6]-[13] を参照．)

問題 2. (n 階導関数の一覧)(1) Table 関数をもちいて，f (x) = x2 sin x の 0 階から 5 階微分までからなるリスト deriv を作成せよ．また，TableForm を適用し一覧表をつくれ．ただし，関数 f[x] の n 階導関数は
D[f[x], {x, n}] で計算できる．
(2) [10] をまねて，deriv の各要素のグラフをまとめて描け．
（描画範囲は各自のセンスで決めてよい．
）
多変数関数の偏微分
[13]
[14]
[15]
[16]

関数 f を再定義：
f のグラフを描く：
x もしくは y で偏微分：：

f[x_, y_] = x^2 + y^2
r = 3; Plot3D[f[x, y], {x, -r, r}, {y, -r, r}]d1 = D[f[x, y], x]d2 = D[f[x, y], y]URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 05-2

担当教員：川平 友規 研究室：A441
[17]
[18]
[19]
[20]
[21]

Email: kawahira@math.nagoya-u.ac.jp

2 階の偏微分，∂ 2 /∂ 2 x：
∂ 2 /∂x∂y ：∂ 2 /∂ 2 y ：f の x 偏導関数 g を定義：
f と g のグラフを描く：（色をつける描画オプション→）：d11 = D[f[x, y], x, x]
d12 = D[f[x, y], x, y]d22 = D[f[x, y], y, y]g[x_,y_] = D[f[x, y], x]Plot3D[{f[x, y], g[x, y]}, {x, -r, r},{y, -r, r}, PlotStyle -> {Blue, Green}])∂2∂2+ 2 f (x, y) = 0 を満問題 3. (ラプラシアン)
∂x2 )
∂y
( 2
∂
1
∂2
∂2
は
たすことを示せ．また，g(x, y, z) = √++g(x, y, z) = 0 を満た∂x2 ∂y 2 ∂z 2
x2 + y 2 + z 2
すことを示せ．
([17][19] を応用して，f (x, y) = ex sin y は問題 4. (極値をもとめる) 関数 f (x, y) = x3 − y 3 − 3x + 12y を考える．

(1) 偏導関数 fx , fy をそれぞれ fx, fy として定義せよ．さらに Solve 関数を用いて，極値を与える点の候補 (sol) を求めよ．また，そこでの f の値を求めよ．
2 で定まる関数を g とおき，sol を代入してその正負を確認せよ．
(2) fxx fyy − fxy(3) Plot3D を用いて，f (x, y) のグラフを描け．ただし，適宜
PlotRange -> {-30, 30}, AxesLabel -> {x, y}などのオプション描画オプションをつけて，極値がどこかわかるように描くこと．
積分

∫
不定積分

[22]
[23]
[24]

sin x dx： int = Integrate[Sin[x], x]グラフを描く： Plot[{Sin[x], int}, {x, 0, Pi}]∫ 2(x2 − 2)dx： Integrate[x^2-2, {x, 0, 2}]
定積分
0

[25]
[26]
[27]

[28]
[29]

積分される部分を図示：
∫ ∞
sin x
広義積分
dx の値：x0
積分される部分を図示：（描画オプションを指示→）：
∫ x
sin t
関数 g(x) =dt を定義：t0
グラフにしてみる：

問題 5. (x の x 乗)

Plot[x^2-2, {x, 0, 2}, Filling -> Axis]
Integrate[Sin[x]/x, {x, 0, Infinity}]
Plot[Sin[x]/x, {x, 0, 10 Pi},
PlotRange -> {-0.5, 1.0}, Filling -> Axis]
g[x_] = Integrate[Sin[t]/t, {t, 0, x}];Plot[{Sin[x]/x, g[x]}, {x, 0, 15 Pi}]x > 0 にたいし，f (x) = xx ，g(x) = f ′ (x)，h(x) =∫ xf (t)dt とさだめる．
0

これらのグラフを 0 ≤ x ≤ 2 の範囲で描け．
（h の計算にはやや時間がかかる．
）

3 次元グラフいろいろ
[30]
[31]
[32]
[33]
[34]

関数 f を再定義：
3 次元のグラフ：
等高線グラフ：
等高線＝陰関数グラフ：
密度プロット：f[x_, y_] = (x^2 + 3 y^2) Exp[1 - x^2 - y^2]
Plot3D[ f[x, y], {x, -2, 2}, {y, -2, 2}]ContourPlot[ f[x, y], {x, -2, 2}, {y, -2, 2}]ContourPlot[ f[x, y] == 1, {x, -2, 2}, {y, -2, 2}]
DensityPlot[ f[x, y], {x, -2, 2}, {y, -2, 2}]

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 05-3

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

問題 6. 問題４の f (x, y) について，ContourPlot および DensityPlot によるグラフを描け．
ContourPlot においては描画範囲やオプション Contours -> 30 などを加え，f (x, y) の極大・極小点および鞍点がわかるように工夫せよ．
課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイ
ル名を 05-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：05-060123456
本文：もし今日の感想・意見・要望などがあれば書いてください．
おねがい：退出時，iMac をログアウトでなくシャットダウンしてください．

研究・ちからだめし
問題 7. (大学入試問題)(1) 2 次関数 f (x) が f (2) = 26, f ′ (2) = 18,∫ 1−1

f (x)dx = 6 をみたすとき，f (x) をもとめよ．

（Hint: f (x) = ax2 + bx + c とおいて a, b, c について Solve）（03 近畿大）∫ x∫ 1(2) 整式 f (x) と実数 C が
f (y)dy +(x + y)2 f (y)dy = x2 + C をみたすとき，この f (x)00と C を求めよ．
（Hint: 等式をよく眺めると，f (x) は 2 次以下であることがわかる．
）（09 京大）(3) x の 3 次関数 f (x) = ax3 + bx2 + cx + d が３つの条件 f (1) = 1，f (−1) = −1,∫ 1(bx2 + cx + d) = 1 をすべて満たしているとする．このような f (x) の中で定積分 I =
−1
∫ 1
2
2
f ′′ (x) dx を最小にするものを求め，そのときの I の値を求めよ．
（Hint: まず 3 条件か−1

ら b, c, d を a の式として Solve．このとき，I は a の関数となる．
）（11 東京大）∫ x∫ 1
1(4) 2 次関数 f (x) と定数 p がf (t)dt+x
f (t)dt− {f (1)−f (−1)} = 4x3 +px2 −10x−4
3
0
−1
を満たす．

(a) 2 次関数 f (x) と定数 p の値を求めよ．
(b) g(x) = (x + 1)f (x) とする．関数 g(x) の極値とそのときの x の値をもとめよ．
（07 同志社大）問題 8. 恒等式をもとめる関数 SolveAlways というのがある．ドキュメントセンターで使い方をしらべて，上の (2)，(4) を解くのに応用せよ．
（もちろん，使わなくても工夫次第で解ける．
）

xy(x2 − y 2 )，h(0, 0) = 0 によって定x2 + y 2
める．この関数は原点で hxy ̸= hyx となる有名な例だが，Mathematica で計算させると hxy = hyxとなることを確認せよ．

問題 9. 関数 h : R2 → R を (x, y) ̸= (0, 0) のとき h(x, y) =URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 06-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

行列と 1 次変換・ParametricPlot
作成日 : November 17, 2011 Version : 1.2

グ ラ フ の ス タ イ ル は「 オ プ ション 」で 調 整 す る ．た と え ば Plot の 場 合 ，Plot のデフォルト値： Plot[Sin[x], {x, -Pi, Pi}]
縦横の縮尺をあわせる： Plot[（同上）, AspectRatio -> Automatic]
x 軸と囲まれる部分を塗る： Plot[（同上）, Filling -> Axis]
グラフの線を緑色で太く： Plot[（同上）, PlotStyle -> {Green, Thick}]オプションはたくさんある（Options[Plot] で一覧が見れる）．使いこなすには，ドキュメントセ
ンター（ヘルプ）を参照しよう．

Tips

:

[a]
[b]
[c]
[d]

◎！注意！◎：問題だけを解くのではなく，[1][2]…等のコマンドも必ず入力し，結果を確かめること．ファイルの提出時もこれらの入力内容は消さずに残しておくこと．

行列とベクトルの演算
[1]

[2]
[3]
[4]
[5]

[6]
[7]
[8]
[9]
[10]
[11]
[12]

ベクトルを定義：
行列を定義：：行列として表示：
ベクトルの内積：ベクトルと行列の積：
行列と行列の積：：行列のべき乗（7 乗）：
行列の転置：
逆行列：
行列式：
トレース：
単位行列・2 次：
3 次：

v1 = {x, y}; v2 = {z, w};
m1 = {{a, b}, {c, d}};m2 = {{3, 2}, {1, 4}};
m1 //MatrixForm もしくは MatrixForm[m1]
v1.v2
m1.v1
pr = m1.m2;
pr // MatrixForm
MatrixPower[m2, 7] // MatrixForm
Transpose[m1] // MatrixForm
Inverse[m1] // MatrixForm
Det[m1]
Tr[m1]
i2 = IdentityMatrix[2]i3 = IdentityMatrix[3](a b問題 1. (ケーリー・ハミルトンの等式) [5][9][10][11] をもちいて，任意の 2 次正方行列 M =c d
2
は等式 M − tr(M )M + det(M )I2 = O をみたすことを示せ．ただし I2 は 2 次の単位行列．

)

固有値・固有ベクトル
[13]
[14]
[15]
[16]
[17]
[18]
[19]

行列を定義：
固有値・固有ベクトルを求める：
対応する固有値・ベクトルは縦に並ぶ：
要素をとりだして…：
定義に当てはめる：
固有多項式：
固有値の確認：m = {{3, 2}, {1, 4}};es = Eigensystem[m]es // MatrixForm
lam = es[[1, 1]]; v = es[[2, 1]];m.v == lam v
p[x_] = Det[m - x i2]
Solve[p[x] == 0, x]

問題 2. (対角化) 上の es をもちいて，q = Transpose[es[[2]]] とおくとこれが m の対角化行列を与える．その理由を考えた上で，Inverse[q].m.q を計算し，確かに対角化されていることを確認せよ．
URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 06-2

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

1 次変換（2 次元）と ParametricPlot
( ) (
)( )
x
3 2
x
2
2
1 次変換 f : R → R を f=で定義する．その作用を可視化しよう．
y
1 4
y
[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]

さっきの行列を確認：
1 次変換を関数として定義：
具体的な値で確認：
円のパラメーター表示：
f による像：
円を表すベクトル値関数：
円とその像を比較：
領域 −1 ≤ x, y ≤ 1：
その像：
比較：
：

m // MatrixForm
f[{x_, y_}] = m.{x, y}v = {-5, 6}; f[v]ParametricPlot[{Cos[t], Sin[t]}, {t, 0, 2 Pi}]
ParametricPlot[f[{Cos[t], Sin[t]}], {t, 0, 2 Pi}]circ[t_] = {Cos[t], Sin[t]};ParametricPlot[{circ[t], f[circ[t]]}, {t, 0, 2 Pi}]
ParametricPlot[{s, t}, {s, -1, 1}, {t, -1, 1}]
ParametricPlot[f[{s, t}], {s, -1, 1}, {t, -1, 1}]
ParametricPlot[{{s, t}, f[{s, t}]},{s, -1, 1}, {t, -1, 1}])
cos mt(0 ≤ t ≤ 2π) で得られる曲線の問題 3. (リサジュー曲線) 整数 m, n を適当に選び，
sin nt
グラフ ParametricPlot で描け．
（整数をいろいろ試し，一番「美しい」と思えるグラフをみつけよ．
）また，その f による像のグラフも描け．
（[25][26] を参考にせよ．
）(

問題 4. (単位円板の像) 関数 disk(r, t) = (r cos t, r sin t) （ただし 0 ≤ r ≤ 1, 0 ≤ t ≤ 2π ）を定
義し，単位円板とその f による像を ParametricPlot で図示せよ．
（[29] を参考にせよ．
）

1 次変換（3 次元）と ParametricPlot3D
  
 
x
1 0 2
x
  
 
1 次変換 g : R3 → R3 を g y  = 0 2 0 y  で定義する．その作用を可視化しよう．
z
2 0 1
z
[30]
[31]
[32]
[33]
[34]
[35]
[36]
[37]
[38]
[39]

行列を定義：：
固有値・ベクトルを確認：
1 次変換を定義：
螺旋 (spiral) を定義：
その 3 次元グラフ：
g による像：比較する：
球面（sphere）を定義：
3 次元グラフとその像：：
描画オプションを指定→：

mm = {{1, 0, 2}, {0, 2, 0}, {2, 0, 1}};mm // MatrixForm
Eigensystem[mm] // MatrixForm
g[{x_, y_, z_}] = mm.{x, y, z}
spi[t_] = {Cos[t], Sin[t], t/10}
ParametricPlot3D[spi[t], {t, 0, 10 Pi}]
ParametricPlot3D[g[spi[t]], {t, 0, 10 Pi}]
ParametricPlot3D[{spi[t], g[spi[t]]}, {t, 0, 10 Pi}]sph[s_, t_] = {Cos[s] Sin[t], Sin[s] Sin[t], Cos[t]};ParametricPlot3D[{sph[s, t], g[sph[s, t]]}, {s, 0, 2 Pi}, {t, 0, Pi},PlotStyle -> {Opacity[0.9, Green], Opacity[0.4, Yellow]}]

問題 5. (トーラス) [38][39] を応用し，以下で定義されるトーラス (torus，円環面) について，そのグラフ，その像のグラフを描け．また，それらを同一空間に描画せよ．
tor[s_, t_] = {Cos[s] (3 + Cos[t]), Sin[s] (3 + Cos[t]), Sin[t]}ただし，s, t の範囲は 0 から 2π までとする．
URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）
担当教員：川平 友規 研究室：A441

演習問題 06-3
Email: kawahira@math.nagoya-u.ac.jp

課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイ
ル名を 06-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：06-060123456
本文：もし今日の感想・意見・要望などがあれば書いてください．
おねがい：退出時，iMac をログアウトでなくシャットダウンしてください．

研究
問題 6. (よりよい可視化) [38][39] の絵は美しいが，1 次変換の作用がわかりづらい．
（たとえば，
この絵からは向きを保つかどうかもわからない．
）そこで，ちょっとした工夫をしてみよう．以下を実行し，1 次変換の作用を確認せよ：[40]
3 枚の板を定義： xy[s_, t_] = {s, t, 0};[41]： yz[s_, t_] = {0, s, t};[42]： zx[s_, t_] = {t, 0, s};[43]
板の絵： ParametricPlot3D[{xy[s, t], yz[s, t], zx[s, t]},：{t, -1, 3}, {s, -1, 3}, PlotStyle -> {Red, Green, Blue}][44] 板の g による像： ParametricPlot3D[{g[xy[s, t]], g[yz[s, t]], g[zx[s, t]]},：{t, -1, 3}, {s, -1, 3}, PlotStyle -> {Red, Green, Blue}]
最後に，[43][44] の絵を同一空間に描け．
問題 7. 2 次元の一次変換 f についても，何らかの工夫をして作用を可視化せよ．
( ) (
)
x
sin(3x+2y)問題 8. 一次変換 f に一見よく似た関数 S : R2 → R2 を S=で定義する．
ysin(x + 4y)
この作用を問題 4 の円板（disk ）関数を用いて可視化せよ．また，半径 r ≤ 0.2 程度の範囲であれば，その絵は f による像とほとんど変わらない．その理由を考えよ．
問題 9. (土星はかっこいい！) 土星（っぽいもの）を描け．

{
問題 10. (大学入試問題) xy 平面上で t を変数とする媒介変数表示

x = 2t + t2で表される曲
y = t + 2t2

線を C とする．

(1) t ̸= −1 のとき，dy
を t の式で表せ．
dxdy1= − を満たす点 A の座標を求めよ．
dx
2

1

 X = √ (2x − y)5(3) 曲線 C 上の点 (x, y) を点 (X, Y ) に移す移動がで表されているとす1 Y = √ (x + 2y)5
る．このとき，Y を X を用いて表せ．
(2) 曲線 C 上で

(4) 曲線 C の概形を xy 平面上にかけ．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/（東京大 06）

名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 07-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

テイラー・フーリエ・Manipulate
作成日 : November 24, 2011 Version : 1.2

Tips : Manipulate 関数の実行結果の中では「Dynamic 変数」というものが使われており，ノートブック中の他の計算と互いに影響しあっている．悪影響を避けるためには，• Manipulate の実行結果は，確認したらすぐに消す（Delete）．
• 変数は頻繁に Clear する．もしくは，できるだけ別の文字を用いるようにする．
◎！注意！◎：問題だけを解くのではなく，[1][2]…等のコマンドも必ず入力し，結果を確かめること．ファイルの提出時もこれらの入力内容は消さずに残しておくこと．

関数の定義（即時割り当て vs. 遅延割り当て）関数を定義する際には，即時割り当てと遅延割り当ての二種類があり，うまく使い分ける必要がある．その違いを，0 と 1 の間の実数をランダムにとる組み込み関数 RandomReal[] を用いて確認してみよう．
[1]
0 と 1 の間のランダムな実数： RandomReal[][2]:= で遅延割り当て： g[x_] := x * RandomReal[]
毎回定義どおり評価される： {g[10], g[10], g[10]}[3][4] 一方，即時割り当ては一旦評価され…： h[x_] = x * RandomReal[][5]
あとは毎回同じ結果となる： {h[10], h[10], h[10]}問題 1. (微分と遅延割り当て) f[x_] = Sin[x]，g[x_]:= D[f[x], x]，h[x_]= D[f[x], x]と定義する．このとき，g[x] は Plot でグラフが描けないが，h[x] はグラフが描けることを確認せよ．また，Evaluate[g[x]] とすれば描けるようになることを確認せよ．

マニピュレート（Manipulate）：動く式
因数分解の Table： Table[Factor[x^n - 1], {n, 1, 10}]
Table を Manipulate に変えてみると…： Manipulate[Factor[x^n - 1], {n, 1, 10}][8]n を整数に限定させる： Manipulate[Factor[x^n - 1], {n, 1, 10, 1}]()
3 2
問題 2. (行列のべき乗) Manipulate を用いて行列 M =のべき乗 M p を p = 1 から 101 4まで整数で変化させよ．
（ちゃんと MatrixForm で見やすくすること．
）
[6]

[7]

マニピュレート（Manipulate）：動くグラフ
[9]
[10]
[11]

[12]
[13]
[14]
[15]
[16]

グラフの Table：
Table を Manipulate に：描画範囲を固定する：:
k を整数に制限：
動く接線を描くぞー！
：
関数・導関数を定義：
接線 (tangent line) を定義：
接線を動かす：::Table[Plot[Sin[k x], {x, -Pi, Pi}], {k, -5, 5}]
Manipulate[Plot[Sin[k x], {x, -Pi, Pi}], {k, -5, 5}]
Manipulate[Plot[(同上), PlotRange -> {-1, 1}],{k, -5, 5}]
Manipulate[（同上）, {k, -5, 5, 1}]
Clear[f, g, x]f[x_] = x^3 - 3 x; g[x_] = D[f[x], x];tanLine[x_, p_] = g[p] (x - p) + f[p]Manipulate[Plot[{f[x], tanLine[x, p]}, {x, -3, 3},
PlotRange -> {-5, 5}], {p, -3, 3}]

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 07-2

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

問題 3. (ガタガタ) [16] で PlotRange -> {-5, 5}] を外すと，グラフがガタガタして情けない感じとなることを確認せよ．
問題 4. (一様収束？) [12] を応用して，fk (x) = xk (k = 1, 2, . . .) がどこで一様収束しているかがわかる「動くグラフ」をつくれ．
問題 5. (接線) [14][15][16] を応用して，好きな関数の「動く接線グラフ」をつくれ．
（他人に見せてそこそこ面白いものがよい．
）

テイラー展開関数を定義：
x=0 での 5 次テイラー展開：
x = Pi/4 での展開：
関数部分を取り出す：
テイラー多項式を遅延的に定義：
即時割り当てだと失敗する：
Evaluate をつけて Plot：
Manipulate で n を変化させる：：：
a も変化させる：
：
：

[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]

[25]

f[x_] = Sin[x];Series[f[x], {x, 0, 5}]Series[f[x], {x, Pi/4, 5}]
Normal[Series[f[x], {x, 0, 5}]]g[x_, a_, n_]:= Normal[Series[f[x], {x, a, n}]]h[x_, a_, n_] = Normal[Series[f[x], {x, a, n}]]Plot[Evaluate[g[x, 0, 5]], {x, - Pi, Pi}]
Manipulate[Plot[Evaluate[{g[x, 0, n], f[x]}],{x, -4 Pi, 4 Pi}, PlotRange -> {-3, 3}],{n, 1, 30, 1}]
Manipulate[Plot[Evaluate[{g[x, a, n], f[x]}],{x, -4 Pi, 4 Pi}, PlotRange -> {-3, 3}],{n, 1, 30, 1}, {a, -2 Pi, 2 Pi}]問題 6. (テイラー展開) [25] を応用して，好きな関数（もちろん Sin 以外の，なんか面白い関数）の「動くテイラー多項式近似」をつくれ．

フーリエ展開区分的に連続な関数 f : R → R および n = 0, 1, 2, . . . にたいし，

1
an =
π

∫ 2π
f (x) cos nx dx,0と定義する．このとき，

f (x) ∼

1
bn =
π

∫ 2π
f (x) sin nx dx,
0

∞
a0 ∑
+
{an cos nx + bn sin nx}
2
n=1

(1.1)

と表し，右辺を f (x) のフーリエ展開 (Fourier expansion) と呼ぶ．
[26]
[27]
[28]
[29]

関数の定義：
関数 f[x] のフーリエ展開：関数として遅延的に定義：
Manipulate で近似の様子を確認：：：

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/f[x_] = x
FourierTrigSeries[f[x], x, 10]g[x_, n_] := FourierTrigSeries[f[x], x, n]Manipulate[Plot[Evaluate[g[x, n]],{x, -2 Pi, 2 Pi}, PlotRange -> {-4, 4}], {n, 1, 100, 1}]
名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 07-3

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

問題 7. [29] を応用して，f (x) = x2 について「動くフーリエ展開」をつくれ．ただし，適当にPlotRange を変えて，大局的な様子がわかるようにすること．
問題 8. (矩形波，Square wave のフーリエ展開) [28] の代わりに次の関数でおきかえて，[29] のように収束の様子を視覚化せよ．(Hint: Sum を用いる．)()
4
sin 3x sin 5x
sin(2n − 1)x
gn (x) :=sin x ++
+ ··· +π35(2n − 1)定理（フーリエ展開の収束）．関数 f : R → R は次の性質を満たすとする：• 周期 2π の関数である．すなわち，f (x) = f (x + 2π)．
• 区間 −π < x < π では f ′ (x) が存在し，連続．
• 極限 f (−π + 0) :=lim

x→−π+0f (x) および f (π − 0) := lim f (x) が存在する．
x→π−0

このとき，連続点 x ̸= (2n + 1)π (n ∈ Z) において∞a0 ∑f (x) =+{an cos nx + bn sin nx}2n=1であり，不連続点 x = (2n + 1)π (n ∈ Z) においては∞
a0 ∑f (−π + 0) + f (π − 0)=+{an cos nx + bn sin nx}22n=1が成り立つ．
（左右からの極限値の平均！）

課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイル名を 07-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：07-060123456
本文：もし今日の感想・意見・要望などがあれば書いてください．
おねがい：退出時，iMac をログアウトでなくシャットダウンしてください．

()
研究cos mt問題 9. (リサジュー曲線) Manipulate と ParamericPlot をもちいて，リサジュー曲線
sin ntが n, m に応じて変化する様子を確認せよ．
問題 10. (等位面の変化) Manipulate と ContourPlot3D をもちいて，陰関数 xyz = k (k ∈ R)で表される 3 次元空間内の曲面が k に応じて変化する様子を確認せよ．
問題 11. (ドーナツを食べる) トーラスを
tor[u_, v_] = {Cos[u] (3 + Cos[v]), Sin[u] (3 + Cos[v]), Sin[v]};で定義する．以下を実行し，トーラスの断面を確認せよ．
Manipulate[ParametricPlot3D[tor[u, v], {u, 0, 2 Pi}, {v, 0, 2 Pi},
PlotRange -> {{-4, 4}, {a, 4}, {-1, 1}}, PlotStyle-> Brown],{a, -4, 4}]

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 08-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

複素級数・複素関数の可視化
作成日 : December 01, 2011 Version : 1.1◎！注意！◎：問題だけを解くのではなく，[1][2]…等のコマンドも必ず入力し，結果を確かめること．ファイルの提出時もこれらの入力内容は消さずに残しておくこと．

今回のテーマは「複素数を視る」ことである．すでに複素関数論を学んだ受講者の皆さんには，それなりの「鑑識眼」があるはず．描くことよりもぜひ，鑑賞するほうに力を注いでほしい．
必要な関数の準備
[1]
[2]
[3]
[4]

複素数から（実部，虚部）のベクトル (vector) へ：
たとえば…：（実部，虚部）のベクトルから複素数 (complex) へ：
たとえば…：

vect[z_] := {Re[z], Im[z]}vect[3 + 4 I]
comp[{x_, y_}] := x + y I
comp[{3, 4}]方程式の虚数解をプロットする
[5]
[6]
[7]
[8]
[9]

方程式の左辺を定義：
方程式を数値的に解く：
ベクトルに直す：
平面上に図示：リストの順番に線で結ぶ：

f[z_] = z^11 - 1;
sol = NSolve[f[z] == 0, z]
v = vect[z] /. sol
ListPlot[v]ListLinePlot[v, PlotMarkers -> Automatic]

問題 1. (より見やすく) [8][9] にオプション AspectRatio -> Automatic を加えて，タテヨコの比率をあわせよ．また，[9] のオプション PlotMarkers -> Automatic を消したらどうなるか確
認せよ．
問題 2.

[5] の方程式を他の高次方程式に変えて，その解を [8][9] のように図示せよ．

Manipulate で Locator を使う
Locator はマウスの動きをそのまま Manipulate のパラメーターとして取り出すことができる
直感的インターフェイスである．
[10]

[11]
[12]
[13]

[14]
[15]

{x,y} 方向の半直線グラフ：：
たとえば…：
Locator で動かす：
0 と z 2 を結ぶ線分：：
たとえば…：
Locator で動かす：

ray[{x_, y_}]:= ParametricPlot[t {x, y}, {t, 0, 5}, PlotRange -> 5]
ray[{-1, 1}]
Manipulate[ ray[v], {{v, {1, 1}}, Locator}]sq[z_] := ParametricPlot[t vect[z^2], {t, 0, 1}, PlotRange -> 5]
sq[1 + 2 I]
Manipulate[ sq[comp[v]], {{v, {2, 1}}, Locator}]

問題 3. (もっと見やすく) [13] のベクトル t vect[z^2] をベクトルのリスト
t {vect[z], vect[z^2]} に変えて，z と z 2 の位置関係が「絶対値は 2 乗，偏角は 2 倍」を満たしていることを確認せよ．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 08-2

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

テイラー級数，およびオイラーの等式 eπi = −1 の可視化
いわゆるオイラーの公式 eiθ = cos θ + i sin θ で θ = π としたものはオイラーの等式とよばれている．円周率 π ，自然対数の底 e，虚数単位 i の豪華メンバーがかくも単純な関係式で結ばれるのだから，不思議なものである．
[16]指数関数 ez を定義： f[z_] := Exp[z];[17] テイラー多項式を関数化： g[z_, n_] := Normal[Series[f[w], {w, 0, n}]]/. w -> z;[18]
z = πi で確認： g[Pi I, 10][19]
Table にしてみる： Table[g[Pi I, k], {k, 0, 10}][20]
ez の近似数列を線分で： seg[z_] := ListLinePlot[
結んだグラフを描く関数：
Table[vect[g[z, k]], {k, 0, 20}],（オプション１）：
PlotRange -> {{-4, 4}, {-4, 4}},
（オプション２）：
AspectRatio -> Automatic,（オプション３）：
PlotMarkers -> Automatic];[21] オイラーの等式を可視化： seg[Pi I][22]z を Locator で動かす： Manipulate[ seg[comp[v]], {{v, {2, 1}}, Locator}]
問題 4. (その他の級数) [16][17][20][22] を応用して，sin z, 1/(1 − z) のテイラー級数が収束もしくは発散する様子を折れ線で「動的に」図示せよ．

複素関数の「グラフ」
[23]
[24]

[25]
[26]

関数・変数を定義：
実部のグラフを描く：：
関数を定義：
z 7→ |g(z)| の 3 次元グラフ：（半径 2 の円板上）：

f[z_] := Sin[z]; z := x + y I
Plot3D[ Re[f[z]], {x, 0, 4 Pi}, {y, 0, 3},
BoxRatios -> Automatic]g[z_] := z^2 - 1; z := r Exp[I t]ParametricPlot3D[{Re[z], Im[z], Abs[g[z]]},{r, 0, 2}, {t, 0, 2 Pi}, PlotRange -> All]問題 5. (複素正弦の虚部と絶対値) [23] を応用して，z 7→ Im (sin z)，z 7→ | sin z| のグラフを描け．
問題 6. (最大値の原理) [26] の出力結果にたいして， g(z) = z 2 − 1 の半径 2 の円板上における最大絶対値は，境界で実現されることを確認せよ．また，この g(z) を ez や sin z に変えても同様の結果が得られることを確認せよ．

複素関数による正方形の像
[27]
[28]
[29]
[30]

[31]
[32]

関数と変数を定義：
正方形の像：
もとの正方形と比較：
w 中心一辺 2r の正方形：とその像を描く関数：：
たとえば…：
Locator で動かす：

f[z_] := z^2; z := x + y I;
ParametricPlot[vect[f[z]], {x, 0, 1}, {y, 0, 1}]
ParametricPlot[{vect[z], vect[f[z]]}, {x, 0, 1}, {y, 0, 1}]
mapf[w_] := ParametricPlot[{vect[z], vect[f[z]]},{x, Re[w] - r, Re[w] + r}, {y, Im[w] - r, Im[w] + r},
PlotRange -> 5]
r = 0.5; mapf[1 + I]
Manipulate[mapf[comp[v]], {{v, {1, 1}}, Locator}]

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）
担当教員：川平 友規 研究室：A441

演習問題 08-3
Email: kawahira@math.nagoya-u.ac.jp

問題 7. (等角性) [32] の結果をよく眺め，もとの正方形の網目の像は，直交する曲線族であることを確認せよ．さらに，もとの正方形の頂点のひとつが原点 z = 0 になるとき，直角が開いて 180
度になること確認せよ．
（必要なら r の値を大きくしてよい．また，mapf に ImageSize -> 500などのオプションを加えると大きな絵が描ける．
）
問題 8. (その他の関数の像) [27][30][32] の関数 f (z) を z 3 , 1/z, ez などに変えて写像の様子を調べよ．これらの場合においても，網目の像の直交性は保たれることを確認せよ．
課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイ
ル名を 08-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：08-060123456
本文：もし今日の感想・要望・質問などがあれば書いてください．できるだけ回答します．
おねがい：退出時，iMac をログアウトでなくシャットダウンしてください．

研究
問題 9. (テイラー級数) [22] の結果に，原点と g[z,0] を結ぶ線分を加えるにはどうすればよいか？
問題 10. (複素関数による円板の像) [32] の正方形を w 中心半径 r の円板に変えよ．また，r も変数として Manipulate できるスライダーを加えよ．
問題 11. (グラフを並べる) [32] の結果において，もとの正方形と像の正方形を別々のグラフに並べて描くにはどうすればよいだろうか？（下の問題がヒントになる．
）
問題 12. (もっと使いやすく) 次を実行せよ．その上で，同等の機能を Locator を用いて実現せよ．
L = 1.3 Pi; f[z_]:= Exp[z];Manipulate[{ ParametricPlot[ {vect[a + I t], vect[t + I b]},{t, -L, L}, PlotRange -> L],
ParametricPlot[ {vect[f[a + I t]], vect[f[t + I b]]},{t, -L, L}, PlotRange -> 5] }, {{a, 1}, -L, L}, {{b, 1}, -L, L}]

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 09-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

確率実験とランダムウォーク
作成日 : December 08, 2011 Version : 1.2◎！注意！◎：問題だけを解くのではなく，[1][2]…等のコマンドも必ず入力し，結果を確かめること．ファイルの提出時もこれらの入力内容は消さずに残しておくこと．

必要な関数の準備
[1]
[2]
[3]
[4]
[5]

1 が 10 個の Table：
Accumulate で足し上げ：
グラフを描く関数：
リストにする：
グラフを重ねあわせる：

tab = Table[1, {10}]Accumulate[tab]gra[n_] := Plot[Sin[n x], {x, 0, Pi}]
gratab = Table[gra[n], {n, 1, 5}]Show[gratab, PlotRange -> All]

乱数を作る関数
[6]
[7]
[8]
[9]
[10]
[11]
[12]

0 か 1 だけの乱数:
たとえば…：|y| ≤ 3 の整数値をとる乱数：
たとえば…：
0 ≤ x ≤ 1 の実数値乱数：
9 ≤ x ≤ 10 の実数値乱数：
リストから等確率で 20 個選ぶ：

x := RandomInteger[];Table[x, {10}]y := RandomInteger[{-3, 3}]
Table[y, {10}]x := RandomReal[]; Table[x, {10}]x := RandomReal[{9, 10}]; Table[x, {10}]RandomChoice[{a, b, c}, 20]

問題 1. (遅延的定義でなきゃだめ) [10] は遅延的定義 x := を用いている．これを即時的定義
x = にすると，結果はどうなるか？
問題 2. (株価っぽいグラフ) 株価 (share/stock)っぽいグラフを描こう．

(1) x := RandomReal[{-1, 1}] とする．dif = Table[x, {100}] にたいし，
shr = Accumulate[dif]] を求めよ．
(2) ListLinePlot[shr] としてグラフを描け．
確率実験（コイン投げ）
コインを 10 回投げて (toss)，表が出た回数だけ金貨（ゴールド）がもらえるとする．このゲー
ムをシミュレーションしてみよう．コインの表を 1，裏を 0 とする．
[13]
[14]
[15]
[16]
[17]
[18]
[19]

コインを 10 回投げる：
たとえば…：
もらえるゴールドをカウント：
たとえば…：
100 人にやらせる：
k ゴールドの人数をカウント：
分布の折れ線グラフ：

toss := Table[RandomInteger[], {10}]tossgold := Total[toss]goldn = 100; trials = Table[gold, {n}]num[k_] := Count[trials, k]
ListLinePlot[ Table[ num[k], {k, 0, 10}] ]問題 3. (ゲームの価格設定) さいころを 5 回投げて，（（出た目の最大値）＋（出た目の最小値））
× 100 円をもらえる有料のゲームを作りたい．1000 人にこのゲームをやらせてみて，損益分岐価
URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 09-2

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

格を自分なりに決定せよ．ただし，リストの総和・最大値・最小値を与える組み込み関数はそれぞれ Total[...], Max[...], Min[...] である．

1 次元ランダムウォーク
数直線の原点に立ち，コインを投げて表なら正の方向へ１移動，裏なら負の方向へ１移動，という操作を繰り返す．その様子をグラフ化してみよう．ただし，コインの表は 1，裏は −1 と考える．
[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]

−1 と 1 をランダムに選ぶ：
n 回選ぶ関数：
たとえば…：
数直線での動きに直す：
たとえば…：
動きをグラフ化：
たとえば…：
500 回を 10 人にやらせる：
グラフを重ねる：

x := RandomChoice[{-1, 1}]p[n_] := Table[x, {n}]p[100]
q[n_] := Accumulate[p[n]]q[100]
r[n_] := ListLinePlot[q[n]]r[10000]
tab = Table[r[500], {10}]Show[tab, PlotRange -> All]問題 4. (グラフの色分け) 重ねたグラフはとても見づらいので，10 人それぞれ別の色になるようにしよう．

(1) 上の [25] で，r[n_, c_] := ListLinePlot[q[n], PlotStyle-> Hue[c/10]] とせよ．
(2) さらに [27] で，tab = Table[r[500, c], {c,1,10}] とし， [28] を実行せよ．

2 次元ランダムウォーク
今度は 2 次元平面の原点からスタートし，「上下左右の方向に，それぞれ確率 1/4 で 1 移動」という操作を繰り返す．その様子をグラフにしてみよう．
[29]
上下左右の方向： dir = {{0, 1}, {0, -1}, {1, 0}, {-1, 0}};[30] 方向を等確率で選ぶ変数： x := RandomChoice[dir][31]
方向を n 回選ぶ関数： p[n_] := Table[x, {n}][32]
たとえば…： p[10][33]p[n] の足し合わせ関数； q[n_] := Accumulate[p[n]]折れ線プロット関数： r[n_, c_] := ListLinePlot[ q[n],[34]（オプション）：
AspectRatio -> Automatic, PlotStyle -> Hue[c/10]];
たとえば…： r[10000, 1]
[35]
[36]
10 人にやらせる： tab = Table[r[5000, c], {c, 1, 10}]グラフを重ねる： Show[tab, PlotRange -> All][37]
問題 5. (偏りのあるランダムウォーク) [29] の dir をうまく修正し，上下に進む確率がそれそれ
1/3 ，左右に進む確率がそれぞれ 1/6 となるようにせよ．そのうえで [36][37] を実行せよ．
課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイ
ル名を 09-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：09-060123456
本文：もし今日の感想・要望・質問などがあれば書いてください．できるだけ回答します．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 09-3

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

おねがい：退出時，iMac をログアウトでなくシャットダウンしてください．

研究：モンテカルロ法 (the Monte Carlo method)
乱数を用いて円周率の近似値を求めることができる：まず，xy 平面上に単位円板を含む正方形{|x|, |y| ≤ 1} をとる．この正方形内にランダムに点をうっていくと，単位円板の中に入っている
点の割合は（円板の面積）÷（正方形の面積）＝ π/4 に近づくであろう．円板に入っているかど
うかは原点との距離の 2 乗で判定できる．また，第一象限だけで実験すればよい．
[38]
[39]
[40]
[41]
[42]
[43]
[44]
[45]

区間 [0, 1] に値をとる乱数：
原点からの距離の 2 乗：
たとえば…：
条件判定関数：
たとえば…：
condi をみたすものを数える：
たとえば…：
円周率を近似する：

x := RandomReal[]; y := RandomReal[];p[n_] := Table[x^2 + y^2, {n}]p[10]condi[x_] := x <= 1;
condi[5]k[n_] := Length[Select[p[n], condi]]
k[100]
n = 10000; N[4 k[n]/n, 5]

研究：3 次元ランダムウォーク

図 1: 3 次元ランダムウォーク

3 次元ランダムウォークも次のようにすれば描くことが出来る．
dir = {{0, 0, 1}, {0, 0, -1}, {1, 0, 0},{-1, 0, 0},{0, -1, 0}, {0, 1, 0}};x := RandomChoice[dir]URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 09-4

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

p[n_] := Table[x, {n}]q[n_] := Accumulate[p[n]]r[n_, c_] := Graphics3D[{Hue[c/20], Line[q[n]]}, AspectRatio -> Automatic]
Show[Table[r[1000, c], {c, 1, 20}], PlotRange -> All]

研究：2 次元ランダムウォークいろいろ問題 6. (六方格子上のランダムウォーク) 正 6 角形のタイルで埋め尽くされた平面を六方格子(hexagonal lattice) とよぶ．[29] の dir と [31] の p[n_] をうまく修正し，六方格子上のランダムウォークを作成せよ．そのうえで [36][37] を実行せよ．(Hint. dir 進んだら次は -dir 進む．)

図 2: 六方格子上のランダムウォーク

問題 7. (動くランダムウォーク) 次を実行すると Manipulate のグラフができる．スライダ横の＋（プラス）を押して，再生ボタンをおすとアニメーションになる．再生速度を遅くして，ランダムウォークの様子を観察せよ．
dir = {{0, 1}, {0, -1}, {1, 0}, {-1, 0}};x := RandomChoice[dir];p[n_] := Table[x, {n}]n = 10000; q = Accumulate[p[n]];qq[k_] := Table[q[[i]], {i, 1, k}]r[k_] := ListLinePlot[qq[k], PlotStyle -> Thick];
total = ListLinePlot[qq[n], PlotStyle -> Red];
Manipulate[Show[{r[k], total}, PlotRange -> All,
AspectRatio -> Automatic, ImageSize -> 500], {k, 1, n}]

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 10-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

ベクトル場と微分方程式
作成日 : December 15, 2011 Version : 1.2◎！注意！◎：問題だけを解くのではなく，[1][2]…等のコマンドも必ず入力し，結果を確かめること．ファイルの提出時もこれらの入力内容は消さずに残しておくこと．

2 次元ベクトル場
ベクトル場 (vector field) v : Rn → Rn とは，Rn の各点 p に「速度ベクトル」 v(p) ∈ Rn を対応付けたものである．風をイメージするのが一番わかりやすい．また，時刻 t ∈ R でパラメーdp(t)ター付けされた曲線 t 7→ p(t) が= v(p(t)) をみたすとき，この曲線をベクトル場 v の流
dt
線 (stream line) もしくは積分曲線 (integral curve) と呼ぶ．
[1]
[2]
[3]
[4]

ベクトル場を描く：
その流線（正弦曲線）：
例その２：
流線は円を描く：

VectorPlot[{1, Cos[x]}, {x, -3, 3}, {y, -3, 3}]
StreamPlot[{1, Cos[x]}, {x, -3, 3}, {y, -3, 3}]
VectorPlot[{-y, x}, {x, -3, 3}, {y, -3, 3}]
StreamPlot[{-y, x}, {x, -3, 3}, {y, -3, 3}]

問題 1. 上と同じ範囲でベクトル場 v = (y, − sin 2x − y) を VectorPlot および StreamPlot せよ．
また，範囲を原点の周りに限定した場合（たとえば |x|, |y| ≤ 0.3），v は u = (y, −2x − y) にき
わめて近いことを確認せよ．それはなぜか？

勾配ベクトル場)∂f ∂f
2 変数関数 f = f (x, y) の勾配ベクトル (gradient vector) とは，ベクトル
,
のことで
∂x ∂y
ある．1 たとえば f が海抜高度を表す関数であるとき，地面にボールを置くと，ボールは勾配ベクトルと「逆向きに」転がっていく．
平面上の関数：
3 次元グラフ：関数の等高線プロット：
勾配ベクトル：
勾配ベクトル場：
その流線：
それらを重ねて表示：
：

[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]

(

f = Cos[x] + y*Sin[x];Plot3D[f, {x, -5, 5}, {y, -5, 5}]
c = ContourPlot[f, {x, -5, 5}, {y, -5, 5}]
gradf = D[f, {{x, y}}]v = VectorPlot[gradf, {x, -5, 5}, {y, -5, 5}]
s = StreamPlot[gradf, {x, -5, 5}, {y, -5, 5}]Show[c, v]
Show[c, s]

問題 2. [5]—[12] を別の関数 g = g(x, y) （自分が面白いと思うもの）に変えて実行せよ．

微分方程式を解く（厳密解と数値解）

Mathematica は（連立）微分方程式の厳密解をDSolve[{方程式，初期条件, ...}, 求める関数名, 変数名]
で与える．
（解けないこともある．
）同様に，数値解は次の形であたえられる：
NDSolve[{方程式，初期条件, ...}, 求める関数名, 変数名とその変域]
1

これを grad f = grad f (x, y) などとあらわす．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 10-2

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

厳密解を求める（正弦曲線）：
解を関数 yy に代入：
値をチェック：
グラフを描く：円の連立微分方程式：：
解をベクトル関数に代入：
グラフ：

[13]
[14]
[15]
[16]
[17]

[18]
[19]

sol = DSolve[{y’[x] == Cos[x], y[0] == 1}, y[x], x]yy[x_] = y[x] /. sol[[1]]yy[Pi]
Plot[yy[x], {x, 0, 2 Pi}]sol2 = DSolve[{x’[t] == -y[t], y’[t] == x[t],x[0] == 1, y[0] == 0}, {x[t], y[t]}, t]p[t_] = {x[t], y[t]} /. sol2[[1]]ParametricPlot[p[t], {t, -Pi, Pi}]

数値解を求めた場合，結果に InterpolatingFunction[...] というのが出てくるが，気にしな
くてもよい．とにかく，厳密解と同様に扱うことができる．
数値解を求める：（ロジスティック曲線）：
解を関数 yy に代入：
値をチェック：
グラフを描く：
線形連立微分方程式：：
解をベクトル関数に代入：
グラフ：

[20]

[21]
[22]
[23]
[24]

[25]
[26]

sol = NDSolve[{y’[x] == y[x]*(1 - y[x]),y[0] == 0.02}, y[x], {x, 0, 10}]yy[x_] = y[x] /. sol[[1]]yy[4]Plot[yy[x], {x, 0, 10}]sol = NDSolve[ {x’[t] == 2 y[t], y’[t] == -2 x[t] - y[t],x[0] == 1, y[0] == 0}, {x[t], y[t]}, {t, 0, 10}];p[t_] = {x[t], y[t]} /. sol[[1]]ParametricPlot[p[t], {t, 0, 10}, PlotRange -> All]

( )
( )
(
)
x′
x
0
2
問題 3. [24] の方程式は=A, A=の形をしている．2 この A を便宜的
y′
y
−2 −1
( )
x
に「定義行列」と呼ぼう．ベクトル場 A
を StreamPlot し，[25] の p[t] がその流線である
y
ことを確認せよ．このような流線はとくに，微分方程式の解曲線もしくは解軌道とよばれる．
問題 4. (初期値を変える) Manipulate の Locator を用いて，[24] の初期値を変えたときの解軌
道の変化がわかるようにしたい．

(1) 次を入力せよ：
str[{qx_, qy_}] := (sol = NDSolve[{x’[t] == 2 y[t], y’[t] == -2 x[t] - y[t],x[0] == qx, y[0] == qy}, {x[t], y[t]}, {t, 0, 10}];p[t_] = {x[t], y[t]} /. sol[[1]];ParametricPlot[p[t], {t, 0, 10}, PlotRange -> 5])
Manipulate[str[q], {{q, {1, 1}}, Locator}]()
02(2) [24] の定義行列 A =を次のものに変えて，(1) と同じことをせよ．また，それ−2 −1ぞれの行列について問題 3 のように StreamPlot し，解軌道の大域的な様子を観察せよ．
(
)
(
)
(
)
(
)
−1 2
4
3
8 −9
2 9
（ア）
（イ）
（ウ）
（エ）
−1 1
−2 11
6 −13
−1 8
2

実は厳密解が存在し，p(t) = p(0) exp(tA) と書ける．ただし，exp(X) = I + X + X 2 /2! + X 3 /3! + · · · .

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）
担当教員：川平 友規 研究室：A441

演習問題 10-3
Email: kawahira@math.nagoya-u.ac.jp

課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイ
ル名を 10-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：10-060123456
本文：もし今日の感想・要望・質問などがあれば書いてください．できるだけ回答します．
おねがい：退出時，iMac をログアウトでなくシャットダウンしてください．

研究：ロトカ-ヴォルテラ方程式
時間 t でのウサギ（被食者）とキツネ（捕食者）の個体数をそれぞれ x = x(t), y = y(t) としたとき，これらの関数は次のロトカ-ヴォルテラ方程式 (Lotka-Volterra Equations) でモデル化で
きる．
x′ = x(α − βy), y ′ = −y(γ − δx)(α, β, γ, δ > 0 は定数)(1) [24][25][26] を応用して，ロトカ-ヴォルテラ方程式を定数 α = γ = 1, β = 0.04, δ = 0.02，
初期値 x(0) = 100, y(0) = 10, 変域 0 ≤ t ≤ 20 として解き，解軌道を xy 平面に図示せよ．
(2) (1) の結果（解軌道が閉曲線のようになる）は何を意味するか考えよ．
(3) (1) の ParametricPlot をただの Plot に変えて結果を比較せよ．
(4) (3) を修正し，Manipulate を用いて δ = 0.02 を 0 ≤ d ≤ 0.1 の範囲で変えられるようにせよ．
（Hint. 問題 4(1) を応用する．この場合 Locator は不要．
）
研究：ローレンツ・アトラクター
次の方程式はローレンツ方程式 (Lorenz equation) と呼ばれる．解曲線が集積する部分はローレ
ンツ・アトラクター (Lorenz attractor) と呼ばれ，極めて複雑な図形となる．
sol = NDSolve[{x’[t] == -10 x[t] + 10 y[t],y’[t] == - x[t] z[t] + 28 x[t] - y[t],z’[t] == x[t] y[t] - (8/3) z[t],x[0] == 1, y[0] == 0, z[0] == 0}, {x[t], y[t], z[t]}, {t, 0, 50}];p[t_] := {x[t], y[t], z[t]} /. sol[[1]];ParametricPlot3D[p[t], {t, 0, 50}, PlotRange -> All]

研究：3 次元勾配ベクトル場

2 次元勾配ベクトル場の 3 次元版．計算量が多いため動作は重い．
[27]
[28]
[29]
[30]
[31]
[32]

3 変数関数：勾配ベクトル：ベクトル場：
等高線：（オプション）：
グラフを重ねる：

f = x^2 + y^2 + z^2;
gradf = D[f, {{x, y, z}}]v = VectorPlot3D[gradf, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]
c = ContourPlot3D[f, {x, -2, 2}, {y, -2, 2}, {z, -2, 2},
Mesh -> None, ContourStyle -> Opacity[0.5, Green]]
Show[v, c]

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 11-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

数列と素数のふるまい
作成日 : December 22, 2011 Version : 1.2

Tips : 下の [1] のように数列や関数の「値」を定義した場合，明示的に上書きしない限り Mathematica
は内部で値を記憶しつづける．悪影響を避けるためには，• 数列・関数・変数を頻繁に Clear する，もしくは• 数列や関数で同じ文字を繰り返し使わない．
◎！注意！◎：問題だけを解くのではなく，[1][2]…等のコマンドも必ず入力し，結果を確かめること．ファイルの提出時もこれらの入力内容は消さずに残しておくこと．

漸化式による数列の定義
[1]
[2]
[3]
[4]

フィボナッチ数列：：
10 項目まで列挙：
折れ線プロット：

Fib[0] = 1; Fib[1] = 1;
Fib[n_] := (Fib[n] = Fib[n - 1] + Fib[n - 2])
tab = Table[Fib[n], {n, 0, 10}]
ListLinePlot[tab]問題 1. (大学入試問題) 数列 {an } を次のように定める．

a1 = 2, a2 = 3, an+2 − an = 4 (n = 1, 2, 3, . . .)このとき，a3 =
40
∑

ak =

， a4 =，a5 =, a6 =

であり，a40 =である．また，である．(Hint. [1][2][3] と同様に 40 項目まで列挙，リストの和は Total[]. )k=1（01 センター本試）

Rsolve で漸化式を微分方程式っぽく解く（いわゆる差分方程式）
[5]
[6]
[7]

2 項間漸化式： sol = RSolve[{a[n + 1] == 2 a[n] + 1, a[0] == 1}, a[n], n]解を数列 aa に代入： aa[n_] := a[n] /. sol[[1]]10 項目までリストにする： tab = Table[aa[n], {n, 0, 10}]問題 2. (再フィボナッチ) [5][6][7] を応用しフィボナッチ数列の漸化式を解き，10 項目までのリストを作れ．
問題 3. RSolve を用いて問題 1 の漸化式を解き，10 項目までのリストを作れ．
関数と Nest による数列の生成
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]

関数 f の 5 回反復：関数を定義：
Nest させる：
数列を定義：
表にする：
表を作る組み込み関数：
たとえば…：
条件関数を決めて…：
その条件で Nest を止める：
たとえば…：

Nest[f, x, 5]
f[x_] := x + 10;
Nest[f, x, 5]a[n_] := Nest[f, 1, n];Table[a[n], {n, 1, 10}]b[n_] := NestList[f, 1, n]b[10]condi[k_] := (k < 200)bb[n_] := NestWhileList[f, n, condi]bb[24]

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 11-2

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

If を用いた場合分け条件に応じて値を決める関数は組み込み関数 If を用いて作ることができる．
If(条件式, True のときの値, False のときの値)
たとえば…： g[x_] := If[x < 0, Sin[x], x]
グラフ： Plot[g[x], {x, -20, 5}]

[18]
[19]

問題 4. (ヘビサイド関数) x < 0 のとき 0，x ≥ 0 のとき 1 となるような関数を作れ．

コラッツの問題与えられた自然数 n にたいし，偶数ならば 2 で割る，奇数ならば 3n + 1 にする，という操作を繰り返すと，かならず · · · 7→ 1 7→ 4 7→ 2 7→ 1 というサイクルに落ち着く，というのがコラッツの予想（the Collatz conjecture, 未解決）である．1
以下では，NestWhileList, If，リストの長さを与える組み込み関数 Length などを用いてこの問題を観察する．
偶数を判定する組み込み関数：
Collatz 関数を If で作成：
たとえば…：「1 でない」ことの条件関数：
操作が終了するまでのリスト：
たとえば…：
グラフ化：
n を動かした表：
TableForm で比較：
操作が終了するのにかかる時間：
n を動かした表：
グラフ化：

[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]

EvenQ[10]col[n_] := If[EvenQ[n], n/2, 3 n + 1];
col[17]
condi[k_] := (k > 1);
colList[n_] := NestWhileList[col, n, condi]
colList[9]ListLinePlot[colList[9]]tab = Table[colList[n], {n, 1, 10}]tab // TableForm
len[n_] := Length[colList[n]]tab2 = Table[len[n], {n, 1, 100}];
ListPlot[tab2, Filling -> Axis, PlotRange -> All]問題 5. 操作が終了するまでの時間が長い n を上のグラフから 1 から 100 までの間で適当に探し，その colList をグラフ化せよ．2
素数の分布
正の実数 x 以下の素数の個数を，整数論では習慣的に π(x) で表す． 関数としての π(x) は階段状の関数であるが，いくつかの近似公式が知られている．たとえば次のものが有名である：xf (x)（ただし f (x) ∼ g(x) :⇐⇒→ 1 (x → ∞) ）
log x
g(x)
∫ x
dt
• チェビシェフの近似：π(x) ∼(この関数は Mathematica の LogIntegral[x])0 log t• ガウスの近似：π(x) ∼10000 番目の素数： Prime[10000]
105000 以下の素数の数： PrimePi[105000]

[32]
[33]
1
2

その他，角谷の問題，3n + 1 予想，シラキューズ予想など，さまざまな呼び名がある．
ちなみに 1 億以下では n = 63728127 が最長らしい．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 11-3

担当教員：川平 友規 研究室：A441
[34]
[35]
[36]
[37]
[38]
[39]
[40]

Email: kawahira@math.nagoya-u.ac.jp

1 万以下の素数のカウント：：
ガウスの近似：：チェビシェフの近似：：グラフを比較：list1 = Table[PrimePi[n], {n, 2, 10000}];
gr1 = ListPlot[list1, PlotStyle -> Red]
list2 = Table[n/Log[n], {n, 2, 10000}];
gr2 = ListPlot[list2, PlotStyle -> Green]
list3 = Table[LogIntegral[n], {n, 2, 10000}];
gr3 = ListPlot[list3, PlotStyle -> Blue]
Show[gr1, gr2, gr3]

問題 6. (近似良し悪し) すでに計算させたリストを用いて，ガウスの近似とチェビシェフの近似の誤差を比較するグラフを描け．(Hint リストはベクトルでもある．
）
課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイ
ル名を 11-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：11-060123456
本文：もし今日の感想・要望・質問などがあれば書いてください．できるだけ回答します．
おねがい：退出時，iMac をログアウトでなくシャットダウンしてください．

研究：双子素数を探す
整数のペア (p, p + 2) がともに素数であるとき，これらを双子素数 (twin prime) とよぶ．双子素数は無限個あると予想されているが，未解決である．ここでは素数のリストから Select 関数をもちいて双子素数を列挙してみよう．
[41]
[42]
[43]
[44]
[45]
[46]
[47]

条件関数：
リストから条件で選択：
双子素数用の条件設定：
素数のペアのリスト：
双子素数を選ぶ：
結構たくさんある：
比率を求める：

condi[n_] := (n > 4);
Select[{1, 4, 6, 7}, condi]
condi[{p1_, p2_}] := (p2 - p1 == 2);
tab = Table[{Prime[n], Prime[n + 1]}, {n, 1, 1000}];
twins = Select[tab, condi];
twins // TableForm
twinsRatio = N[Length[twins]/Length[tab]]問題 7. (4k + 1 vs 4k + 3) 2 以外の素数は一般に 4k + 1 もしくは 4k + 3 の形をしている．上の双子素数を見つける操作を参考に，4k + 1 および 4k + 3 の形の素数をそれぞれ取り出してリストにせよ．また，それぞれの比率を求めよ．ただし，n を 4 で割ったあまりは組み込み関数 Mod[n, 4]で与えられる．
問題 8. (三つ子素数 prime triplet) 整数の三つ組 (p, p + 2, p + 6) または (p, p + 4, p + 6) がす
べて素数からなるとき，これを三つ子素数とよぶ．三つ子素数を列挙せよ．三つ子素数も無限個あると予想されているが，未解決である．ただし，「x = 0 かつ y = 0」は (x == 0) && (y == 0)，「x = 0 または y = 0」は (x == 0)||(y == 0) のように表す．

研究

ϕn+1
問題 9. (偽フィボナッチ) 数列 bn = √ (n ≥ 0) はフィボナッチ数列を極めてよく近似するこ5√とを確認せよ．ただし，ϕ = ( 5 + 1)/2 は黄金比であり，Mathematica では GoldenRatio と表される．
URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 11-4

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

問題 10. (ロジスティック・カオス) Nest を用いて数列のカオス的なふるまいを観察しよう．以下を入力すると，関数 fc (x) := cx(1 − x) により 1/2 7→ fc (1/2) 7→ fc (fc (1/2)) 7→ · · · で与えられる数列が c に応じてカオス的に変化することが観察できる．
f[c_][x_] := c x (1 - x);lis[n_, c_] := NestList[f[c], 0.5, n]
gra[n_, c_] := ListLinePlot[lis[n, c], PlotRange -> {{0, n}, {-2, 2}}]
Manipulate[gra[20, c], {c, -5, 5}]
問題 11. (バンプ関数) 多様体の微積分で重要な C ∞ 関数である，バンプ関数 (bump function)のグラフを描こう．
{
0(x ≤ 0)を定義し，グラフを描け．
(1) If をもちいて，a(x) :=2−1/xe(x > 0)(2) b(x) :=a(x)のグラフを描け．
a(x) + a(1 − x)(3) 適当に ϵ > 0 をとり，c(x) := b(x/ϵ + 2)b(−x/ϵ + 2) のグラフをかけ．
問題 12. (続・素数の分布) 歴史上有名な π(x) の近似式として，x
log x − 1.08366
)
∞ (
∞
∑
∑
(log x)n1• リーマンの近似：π(x) ∼ 1 +（ただし ζ(s) =）n · n! · ζ(n + 1)
ns
• ルジャンドルの近似：π(x) ∼

n=1

n=1

が知られている．これらの近似式と π(x) をグラフ上で比較せよ．また，ガウス，チェビシェフのものとも比較せよ．ただし，リーマンの近似式は組み込み関数 RiemannR[x] で与えられる．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 12-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

セル・オートマトン
作成日 : January 12, 2012 Version : 1.1

ArrayPlot による行列の表現
ArrayPlot は行列の各要素を正方形のセル（もしくはピクセル）で表現する．ただし，ここでいう「行列」とは，「リストのリスト」のこと．
[1] 0 と 1 によるランダムな 10 × 20 行列： m = Table[RandomInteger[], {10}, {20}];
TableForm で表現： TableForm[m][2][3]ArrayPlot で表現： ArrayPlot[m][4]0 と 1 の間のランダムな実数の場合： m = Table[RandomReal[], {10}, {20}];[5]： TableForm[m][6]： ArrayPlot[m]問題 1. (カラーにする) 0 と 1 の間の実数 x にたいし，Hue[x] とすれば色相環1 に基づいた「色」が指定できる．たとえば Hue[0] は Red，Hue[1/3] は Green，Hue[2/3] は Blue といった具合である．
いま行列を m = Table[Hue[RandomReal[]], {10}, {20}]; と定めて，これを TableForm お
よび ArrayPlot で表せ．また，行列のサイズを適当に変化させよ．

パスカルの三角形n!
を並べたものk!(n − k)!
である．これは (x + 1)n を展開したときの xk の係数に等しく，さらに関係式
パスカルの三角形 (Pascal’s triangle) とは，いわゆる 2 項係数 n Ck =n+1 Ck= n Ck−1 + n Ck

が成り立つのであった．これらの数を整数で割ったあまりがどのように分布するか，ArrayPlot
で観察してみよう．ただし，p を q で割ったあまりは組み込み関数 Mod[p,q] で与えられる．

2 項係数を定義：一覧表にする：
TableForm で表現：
n=100 のとき，mod 2 の表：ArrayPlot する：

[7]
[8]
[9]
[10]
[11]

c[n_, k_] := n!/(k!*(n - k)!)pas = Table[c[n, k], {n, 1, 10}, {k, 0, n}];TableForm[pas]pas2 = Table[Mod[c[n, k], 2], {n, 1, 100}, {k, 0, n}];ArrayPlot[pas2]

問題 2. (カラーにする) pascal[j_] = Table[Hue[Mod[c[n, k], j]/j], {n, 1, 100}, {k, 0, n}];とおき，Manipulate を用いて，pascal[j] を ArrayPlot したものが j が 2 から 100 まで変化できるようにせよ．
（j は整数値のみを取ることに注意！）

If による条件分岐
If で不連続関数を場合分けで定義：
そのグラフ：
If を繰り返し，さらに場合分け：
：

[12]
[13]
[14]
[15]
1

f[x_] := If[x < 0, -1, 1];
Plot[f[x], {x, -2, 2}]g[x_] := If[x < 0, -1, If[x < 1, 1/2, 1]];
Plot[g[x], {x, -2, 2}]

赤〜緑〜青〜赤という色の変化が円状に配置されたもの．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）
担当教員：川平 友規 研究室：A441

演習問題 12-2
Email: kawahira@math.nagoya-u.ac.jp

問題 3. (周期性をもった数列の作成) If と Mod[n,3] を用いて，π, e, i, π, e, i, π, . . . と続く
数列 a[n] を a[n_] := If[...] の形で定義し，100 項目までのリストを作成せよ．
セル・オートマトンを自作する
周期的境界をもつセル・オートマトン (cellular automaton) を作成しよう．
（詳細は講義で説明する．
）[16] 初期状態の作成（100 個の 0）： n = 101; c0 = Table[0, {n - 1}];
[17]
50 番目に 1 を挿入： c = Insert[c0, 1, (n - 1)/2];
ArrayPlot で表現： ArrayPlot[{c}]
[18]
[19]
ルール 150 を定める関数： f[x_, y_, z_] := Mod[x + y + z, 2];[20]
i 番目のセルに f で定まる： rule[c_, i_] :=ルールを適用する関数：
If[i == 1, f[ c[[n]], c[[1]], c[[2]] ],：If[i == n, f[ c[[n-1]], c[[n]], c[[1]] ],：f[ c[[i-1]], c[[i]], c[[i+1]] ](ふたつの If の閉じ括弧→)：]]次のセル状態を作る関数： next[c_] := Table[rule[c, i], {i, 1, n}];
[21]
[22]
もとの c に適用： c = next[c]; ArrayPlot[{c}][23]
もう一回： c = next[c]; ArrayPlot[{c}][24]さらにもう一回： c = next[c]; ArrayPlot[{c}][25]
その続きをタテに並べる： ArrayPlot[NestList[next, c, 200]]問題 4. (その他のルール) 以上を応用して，次のように変化させよ．

(1) [18] の初期状態から，[20] を変形させた関数 f[x_, y_, z_] := Mod[x + y + z + y z, 2];によるルールを適用した場合，0 ステップ目から 200 ステップ目までの状態を [26] のように
ArrayPlot せよ．
（これはルール 30 となる．
）
(2) (1) の初期状態を 0 と 1 からなるランダムな列に変えた場合はどうか？
(3) ほかにもルールを自作して，どのようなパターンが得られるか確認せよ．
（状態数を 0, 1, 2などに増やすのも面白い．Hue を使ってカラーにしても面白い．この場合，一旦 200 ステッ
プ目までの NestList を作ってから，各要素が 0 と 1 の間の値になるように調整し，さらに各要素の Hue からなる行列を Table で作るとよい．
）

組み込み関数によるセル・オートマトン

Mathematica の開発者 S. Wolfram はもともと，セル・オートマトン理論のエキスパートである．もちろん Mathematica の中にも，セル・オートマトンを生成するための組み込み関数が準備されている．
[26]
[27]
[28]

初期状態を作成： c0 := Table[0, {100}]; c := Insert[c0, 1, 50];c にルール 30 を 200 回適用： b := CellularAutomaton[30, c, 200];
ArrayPlot で表現： ArrayPlot[b]問題 5. (いろんなルール) Manipulate をもちいて，初期状態 c にルール 0 から 255 までを適用
したセル・オートマトンが図示できるようにせよ．
（ルールの番号は整数だけをわたるようにする
こと．
）
URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 12-3

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイ
ル名を 12-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：12-060123456
本文：もし今日の感想・要望・質問などがあれば書いてください．できるだけ回答します．
おねがい：退出時，iMac をログアウトでなくシャットダウンしてください．

図 1: イモ貝の殻（出典：Wikipedia）

研究
問題 6. (アニメーション) Manipulate を用いて，ルール 30 によるセル・オートマトンのアニメーションを作成せよ．
（これにはふたつの方法が考えられる．一列のセルが時間に応じて変化するものと，列が下に伸びていって時間発展の様子がわかるように変化するもの．いずれでもよい．
）

問題 7. (いろんな境界条件) 上であつかった n 個のセルからなるセル・オートマトンでは，境界がなく，セルは周期的（すなわち円環状にセルが並んでいる）と考えた．これは列の両端に仮想的なセル c[[0]] = c[[n]] ，c[[n + 1]] = c[[1]] を考えていることに対応する．そのほか，
境界条件として次のものがある．

• 固定的： 仮想的に c[[0]] と c[[n + 1]] を定数として設定する．
• 反射的： 仮想的に c[[0]] = c[[2]]，c[[n + 1]] = c[[n-1]] として設定する．
これらの場合，セル・オートマトンがどのように変化するか確かめよ．
（n が大きすぎると，変化がわかりづらいかもしれない．
）

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 13-1

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

フラクタル
作成日 : January 19, 2012 Version : 1.2

Nest 以外のループ処理（Do, For, While）
関数や漸化式の反復は Nest を用いて実現できた．それ以外にも，Mathematica には一般的なプログラミング言語で使われるようなループ処理（反復処理）用の組み込み関数が備わっている．
以下はすべて，1 から 5 までをタテに並べて表示するためのプログラムである：
[1]
[2]
[3]

Do の使い方： Do[Print[i], {i, 1, 5}]
While の使い方： i = 1; While[i < 6, Print[i]; i = i + 1]
For の使い方： For[i = 1, i < 6, i = i + 1, Print[i]]これらは次のように解釈する：• Do は，「Print[i] を i = 1 から i = 5 まで（1 刻みで）動かして実行せよ」．i の刻み幅を 0.5 にしたければ，Do[Print[i], {i, 1, 5, 0.5}] のように指定できる．
• While は，「まず i = 1 とせよ．つぎに条件 i < 6 が満たされている限り，Print[i] とi = i + 1 を実行せよ」• For は，
「まず i = 1 とせよ．つぎに条件 i < 6 が満たされている限り，i = i + 1 とPrint[i] を交互に実行せよ」1
問題 1. (プログラムを読む) 『x = 1; Do[x = x * i, {i, 1, 5}]; x』と『x = 1; Do[x = x * i, {i, 5, 1, -1}]; x』はともに 5! を計算していることを納得せよ．また，同等のプログラムを While を用いて書け．

リストをほどく，置換する
フタクタル図形を描くにあたって，必要なリストの操作を理解しておこう．
リストのほどき方（平坦化）
[4]
[5]
[6]

3 重のリスト： list = { {{a, a}, {b, b}}, {{c, c}, {d, d}} };全レベル平坦化： Flatten[list]
1 レベル平坦化： Flatten[list, 1]

リストの和集合[7]{a , b} と {c , d} の和集合： Join[{a , b}, {c , d}]

リストの要素の置換[8]

3 番目の項を xxx に置換： ReplacePart[{a, b, c, d, e}, 3 -> xxx]

リストの各項に関数を施す[9]

リストの各要素に関数 f を適用： Map[f, {a, b, c, d}]1

C 言語など他のプログラミング言語と比べて，For の使い方が若干ややこしい．プログラムの読みやすさを考えると，Mathematica では使わないほうがいいかもしれない．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 13-2

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

コッホ曲線
複素平面上で，0 と 1 を結ぶ線分からスタートし，次のような操作を限りなく続けて得られる「極限」の図形がコッホ曲線 (Koch curve) である．
（「極限」の数学的な意味，詳しい構成方法については授業で解説する．
）この図形を（近似的に）作図しよう．

関数の準備． いくつか必要な関数を定義しておく．
[10]
[11]
[12]
[13]
[14]

[15]

複素数→ベクトル：
たとえば：複素数のリスト→ベクトルのリスト：
たとえば：複素数のリストを線分で結ぶ：
：
：
：
たとえば：

vect[z_] = {Re[z], Im[z]};vect[1 + I]
vectorize[list_] := Map[vect, list]
vectorize[{-1, I, -I, 1}]complexLLP[list_] :=ListLinePlot[vectorize[list],
AspectRatio -> Automatic,
Axes -> None, ImageSize -> 600];
complexLLP[{-1, I, -I, 1}]

問題 2. (Map の効用) vectorize2[list_] := vect[list] としたとき，vectorize2[{-1, I, -I, 1}]
はどうなるか？（これでは ListLinePlot を適用しても思うような絵にならない！）

コッホ曲線を描く．
[16]
[17]

[18]
[19]

定数を決める：
リストの置換用関数：：
最初の線分：
次の線分（ループ）：：：：：(図示する)：

s = Exp[Pi/6 I]/Sqrt[3.0] ;rep[a_, b_] :={a + (b - a)/3, a + s*(b - a), a + (b - a)*(2/3), b};koch = {0.0, 1.0};
koch = (ko = koch;Do[ko = ReplacePart[ko, i -> rep[ ko[[i-1]], ko[[i]] ]],{i, Length[koch], 2, -1}];
Flatten[ko]);complexLLP[koch]

問題 3. (手動で反復) [19] を繰り返し実行して（再入力は不要だよ，念のため），コッホ曲線に近づいていく様子を観察せよ．
問題 4.

[19] を While を用いて書き直せ．
（時間が足りなそうなら，省略して次にすすむこと．
）

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）

演習問題 13-3

担当教員：川平 友規 研究室：A441

Email: kawahira@math.nagoya-u.ac.jp

IFS によるコッホ曲線
eπi/6
定数 s = √ を用いて関数 Ti : C → C (i = 1, 2) を
3
T1 (z) = −s(z − 1),T2 (z) = −sz + 1

で定めるとき，コッホ曲線は• K0 = (0, 1, s を頂点とする三角形)• Kn+1 = T1 (Kn ) ∪ T2 (Kn )として得られる図形の列 {Kn }n≥0 の「極限」とも考えられる．

[20]
[21]
[22]
[23]
[24]

さっきの定数：関数その１：
関数その２：K0 ：
Kn から Kn+1 へ：（複素数→ベクトル）：
図形の描画：

s = Exp[Pi/6 I]/Sqrt[3.0] ;T1[z_] := -s (z - 1);
T2[z_] := -Conjugate[s] z + 1;
ver = {{0, 1, s}}ver = Join[Map[T1, ver], Map[T2, ver]];triangles = Map[vectorize, ver];
Graphics[Polygon[triangles], ImageSize -> 600]問題 5. (手動で反復) [24] を繰り返し実行して，コッホ曲線に近づいていく様子を観察せよ．
問題 6. (シェルピンスキー・ガスケット)

S1 (z) = z/2,

定数 α = eπi/3 を用いて関数 Si : C → C を

1
S2 (z) = (z − 1) + 1,
2

1
S3 (z) = (z − α) + α2で定めるとき，• K0 = (0, 1, α を頂点とする三角形)• Kn+1 = S1 (Kn ) ∪ S2 (Kn ) ∪ S3 (Kn )として得られる図形の列 {Kn }n≥0 の「極限」はシェルピンスキー・ガスケット（Sierpinski’s gasket）と呼ばれる．上の [20]—[24] を応用して，K1 , K2 , . . . を描画せよ．

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）
担当教員：川平 友規 研究室：A441

演習問題 13-4
Email: kawahira@math.nagoya-u.ac.jp

課題の提出：以上の問題を解いたノートブック（.nb）ファイルをパッケージ形式（.m）として保存し，ファイ
ル名を 13-060123456.m のようにしてメールで提出する．
あて先：comp0-2011@math.nagoya-u.ac.jp
件名（半角で！）：13-060123456
本文：もし今日の感想・要望・質問などがあれば書いてください．できるだけ回答します．
おねがい：退出時，iMac をログアウトでなくシャットダウンしてください．

研究
問題 7. (ドラゴン) 次のようにして得られる曲線の「極限」はハイウェイ・ドラゴン (Heighway’s
dragon) と呼ばれる．これを描画せよ．

問題 8. (カーペット) 正方形から中央の正方形を次のように除いていって得られる「極限」はシェ
ルピンスキー・カーペット (Sierpinski’s carpet) と呼ばれる．IFS を用いてこれを描画せよ．

問題 9. (ペンタ君？) 正五角形からスタートし，次のような操作を繰り返して得られる図形の「極限」はペンタクン (pentakun) と呼ばれる．IFS を用いてこれを描画せよ．(Hint: 1 の 5 乗根を頂√点にもつ 5 角形からスタート．これらの頂点を中心に複素平面を縮小率 (3 − 5)/2 で縮小する写像を用いる．)

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科

計算数学基礎（2011 後期）
担当教員：川平 友規 研究室：A441

演習問題 13-5
Email: kawahira@math.nagoya-u.ac.jp

問題 10. (カラーでペンタ君) 上のようにランダムなカラーでペンタクンを描け．(Hint: 色のついた多角形は Graphics[{Red, Polygon[..]}] のように指定する．)
参考． 今回紹介したような図形はフラクタル (fractal) と呼ばれる．これらの共通する性質は自己相似性 (self-similarity) である．このような図形には整数でない次元を定義することができる．たとえば…
4• コッホ曲線： log
log 3 ≈ 1.262 次元
3• シェルピンスキー・ガスケット： log
log 2 ≈ 1.585 次元• ハイウェイ・ドラゴン：2 次元（平面を部分的に埋め尽くす）
8• シェルピンスキー・カーペット： log
log 3 ≈ 1.893 次元√ 5• ペンタクン： log(3−log≈ 1.672 次元5)−log 2

URL: http://www.math.nagoya-u.ac.jp/˜kawahira/名古屋大学・理学部・数理学科