数学あいかわらず苦手だけど、数学得意な人達にとってはプログラミング言語と同様の便利なツール、って感覚が分かってきた。 https://t.co/lgHCS0IVTY pic.twitter.com/cvZ1q8WEsv
— Kazunori Sato (@kazunori_279) 2016年10月4日
参考にするためまずはメモ.
どういうコンテンツ作っていこう?
とりあえずは Sympy 頑張ろう.
数学あいかわらず苦手だけど、数学得意な人達にとってはプログラミング言語と同様の便利なツール、って感覚が分かってきた。 https://t.co/lgHCS0IVTY pic.twitter.com/cvZ1q8WEsv
— Kazunori Sato (@kazunori_279) 2016年10月4日
参考にするためまずはメモ.
どういうコンテンツ作っていこう?
とりあえずは Sympy 頑張ろう.
こちらに PDF があります.
サイトでは見づらい方, コンテンツを手元に置いておきたい方は
ダウンロードしてご覧ください.
最初に偏微分方程式のシミュレーション・アニメーションを見てもらいました.
物理, 経済, 生物で出てくる微分方程式を紹介し,
それを差分近似でシミュレーションしました.
このシミュレーションの中で出てくる数学にベクトル,
関数, 数列, 漸化式, そしてもちろん微分があります.
それらを順に説明してきました.
高校で学ぶ関数の大事な具体例として指数関数,
三角関数や対数関数があります.
指数関数と三角関数について微分方程式を解く最初の段階で紹介しました.
中学から出てくる一次方程式については,
シミュレーションで実際に解いたのが一次方程式ですし,
微分の基礎に横たわっているのが一次近似という発想だ,
ということも紹介しました.
説明不足なところも多いです.
しかし細かいところばかりにこだわって大局を見渡せないのが問題だ,
というのがこの講座を作った動機です.
まずはこの大きな流れ,
社会を支え社会に生きる数学の姿を感じてもらうことに集中しました.
意識するところは大きくわけて 2 つあるでしょう.
これはもちろんあなたが次のどのタイプに入るかによります.
どんな興味があるかによっていったん道がわかれます.
その道はあとで合流するかもしれないし,
合流せずとも細く長くずっと続いていくかもしれません.
タイプごとに指針や参考書を紹介します.
都合のいい部分をうまく使ってみてください.
これからいろいろな本やコンテンツを紹介していきます.
しかしその前にこの講座に参加したあなたに必ず読んでほしい文章があります.
もともと山形大学の数学エッセイコンテストで入賞していた作品です.
理学部や数学科の HP 改訂でどこにあるのかわからなくなってしまいました.
もちろん著作権などの微妙な問題もあるのですが,
埋もれさせるにはあまりにも惜しい文章です.
そしてこれは何をさしおいても真っ先に読んでください.
「数学とこんな付き合い方もできるのか」とはっとさせられる,
穏やかな筆致のとても素敵な文章です.
6 ページの短い文章なのでちょっとした隙間時間でも読み切れます.
簡単にあらすじを言うと,
中学の数学で挫折した女性がお子さんと一緒に算数から勉強をはじめ,
数学の呪縛から解き放たれたという文章です.
勉強の仕方もとても示唆的です.
小学校 1 年からお子さんとその日習ったことを一緒に勉強していくスタイルです.
要は自分で「今日はこのくらい」としたのではなく,
お子さんの授業の進度という他人のペースに合わせて,
9 年かけて中学数学まで勉強を続けた事例です.
子どもからすれば全科目で毎回新しいことをたくさん勉強するわけで,
とんでもないハードワークです.
しかし算数だけ取り出すなら,
それも大人の視点から見るのなら,
十分小さなボリュームに小分けされていると思っていいでしょう.
今回の通信講座のスタイルはもちろんこれを参考にしています.
大人の勉強スタイルはいくつかあります.
何をどうしたいかによって使い分ける必要があります.
そしてあなたの目標によっては短期決戦スタイルの勉強が不向きな場合があります.
この植村さんの事例を単なるストーリーとして消化するのではなく,
ここで語られた無理のない勉強スタイルをどう自分の生活に組み込むか,
そこまで考えながら読むようにしてください.
「こんな本がいいですよ」とか,
「こんなサービスがいいですよ」とかいう具体的な話をする前に,
『たかが数学, されど数学』での大事な教訓をまとめます.
端的に言えば 1 人で勉強するのはつらいから誰かを巻き込みましょう,
ということです.
1 人で勉強できるなら最初から中高数学で困ってないはずです.
本の選び方がよくないとかいう話よりも,
こちらの方が本質的です.
そもそも本当に数学を勉強したいんですか?というのもけっこう微妙な問題です.
この辺を掘っていきましょう.
環境の重要性は勉強・学業という視点からはあまり強調されません.
しかしとても大事なことです.
わかりやすいのは部活やスポーツでしょうか?
「ちゃんとやるならいい環境で」というのは,
いたってふつうの考えとして浸透しているように思います.
クラブチームや強豪校からプロへ,
こうした流れの中で他人と切磋琢磨する環境が大事だ,
とはよく言われることでしょう.
これは勉強する上でも同じです.
私は中高はふつうの公立校でした.
しかし大学に入っていわゆる進学校と呼ばれる学校に通っていた友人達の話を聞いたら,
環境というか世界が違うことに気付きました.
勉強するのは当たり前だし,
いろいろな情報が入ってくるし,
何より楽しそうに勉強している人が周りに多いから自分も勉強していて楽しくなる,
そういう常識・感覚レベルで全く違う世界に生きていました.
ちなみに私, 高校で文系理系を選択するとき,
理系の人ってみんな数学好きだと思っていたんですね.
進級したら全然違うのでめちゃくちゃ驚いて,
「こんなのは自分だけなのか」と思っていました.
大学に入ってようやく「ああ,
同類はこんなにたくさんいるじゃないか」とわかりました.
環境が違うと本当に世界や感覚, 常識が全然違います.
意図的かどうかはさておき,
植村さんはすさまじい環境を整備し,
尋常ならざる目標を立てています.
目標は「小学校に上がった子どものペースに合わせて 9 年かけてやり直す」です.
3 日坊主という言葉があるくらい,
行動が長続きしないのがふつうです.
その中で 1 年どころか 9 年がかりでやり通すというのは尋常なことではありません.
そして環境.
お子さんの学習状況というペースメーカーがあります.
毎日たくさんの科目で新しいことを習い続ける子どもからすればしんどくても,
算数に的を絞った大人からすればまあゆっくりです.
無理が全くないとも言えます.
しかも次のような記述もあります.
知りたがり屋のオカンがすぐ頭を突っ込んでくるので面白がって今日習ったところを得意になって教えてくれた.
母親が面白がって聞いてくる.
子どもたちも面白がって楽しく勉強する.
わからないところは得意になって教えてくれる.
お子さんは無料の教師になってくれてさえいます.
はっきり書いてないですが,
植村さん, 自分がわからないところをわかるように教えてくれたら,
自分のお子さんをすごい褒めたと思うんですね.
そういういいフィードバックのループが回る環境を作ったわけです.
勉強が続かないわけがないですね.
これが, この環境構築が本当に大事.
他人の力を借りるというのは,
簡単に言うとそういう環境を作ることです.
「何か勉強しないとまずいよな」とか,
「私も何か勉強してみたい!」,
そう思える環境を作ってそこに所属することです.
要は 1 人でやっても続かない現実とどう向き合うか,
そういう問題です.
植村さんのようにお子さんがいるならそれを刺激にするのもいいでしょう.
お子さんは教師にもペースメーカーにもなってくれます.
お金もかかりませんしね.
もちろんこんな都合のいい環境が自然に作れる人はいないでしょう.
お子さんが「勉強嫌い」ならそもそもお子さんの力も借りられないでしょうし.
そうなると他で環境を作る必要があります.
最近あまり聞きませんが,
特にあなたが都会にお住まいなら,
朝活に顔を出して毎日朝 30 分は勉強する,
そういう環境を作ったり使ったりしてもいいでしょう.
みんな何か目的に向かって行動している環境に身を置けば,
自然にあなた自身も行動できるようになれますよ.
ただ, これはこれで割と強い意志が必要です.
朝活の場に行き続けないといけないですから.
他人を本気で巻き込もうと思うなら,
積極的に話しかけたり勉強会の場を作ったりするのも必要ですしね.
そういう観点からすると,
1 番手っ取り早いのは児童・生徒向けの塾などに参加することです.
時間が合わないなら進研ゼミなどの通信教育もあります.
ちなみに私の知人の女性で「彼氏からもらった数学ガールをちゃんと読んでみたくて」と言って,
公文式に通って教材をもらい,
自分でコツコツ進めている人がいます.
「因数分解って楽しいですよね!」と本当に楽しそうに言っていて,
「ああ, いいな」と素直に思いました.
はっきり言えば 1 人でやっていても心が折れます.
だから他人を巻き込もう, そういう話です.
こういうことを言うと嫌がる人も多いのですが,
お金を払って環境を作るのはお勧めですよ.
それは自分の覚悟の証でもあります.
逆にそうまでしてもやれないなら,
それはあなたにとって算数や数学なんて必要ないからです.
継続的にお金を払ってまでやりたいか,
これは判断基準としてかなり使えます.
これ, 本当に真剣な話なんですが,
例えば数学できなくて困ったことありますか?
学生の頃に試験で困った以外の経験ありますか?
どうしても数学したいというのに何で自分 1 人で続けられないんですか?
もしも好きなことなら意志の力なんて言う余地すらなく勝手に続けますよね?
その辺までちゃんと自分の心と相談してみてください.
環境を作るのが大事とは言いました.
しかしなかなかその環境を作るのも難しいのが現実です.
要は独学せざるを得ないこともよくあります.
いま本についても調べているところです.
しかしある程度の薄さで中高全範囲をカバーしつつ,
発展的な話題も程々に扱っている都合のいい本はあまり見かけません.
発見次第適宜紹介したいとは思っています.
くり返しになりますが,
本の選定基準を改めてまとめます.
この 2 点を重視するのがお勧めです.
少し違う視点からのお勧めがあるので,
それも紹介しましょう.
それはストーリーに数学が埋め込まれた本を読むことです.
本当に数学もきちんと勉強できるという点からすると次の数学ガールがお勧めです.
ここでは一冊しか紹介していませんが,
ちょうど中高の数学に対応するレベルの内容が
「秘密のノート」シリーズとしてシリーズ化されています.
1 冊 1 冊は 200 ページ程度あるし,
全部買うとそれなりの値段にもなります.
しかしそれだけのお金をかける価値のある本・シリーズだとも思います.
次の本はいわば本編, 大学の数学にも足を踏み入れた内容です.
これもかなり発展的な話題を丁寧に追いかけています.
特定の分野を詳しく追いかけるというより,
ある問題を考えると自然に関係してくるいろいろな数学の世界をのんびり垣間見ていくという感じです.
こちらもシリーズで既に 5 冊以上出ているはずで,
全部買うとそれなりの値段になります.
しかしこれもそれだけの価値はあります.
あなたはもっと気楽に数学を楽しみたいと思っているかもしれません.
息抜きに数学をネタにした小説を読んでみたい,
そんな要望もあるでしょう.
その手の本やコンテンツは最後にお伝えします.
ご興味があればそちらも見てみてください.
大人向けの数学に関するリアルの教室を 2 つ知っているので,
そちらもご案内しておきます.
どちらも東京にあります.
和は大阪にも展開しているようです.
これは私の活動と全く関係ないサービスです.
評判は悪くないみたいなのでとりあえずご紹介, という感じです.
両方とも有料でそれなりに値段も張ります.
リアルでじっくり質問したいといったご要望がある方にはいいでしょう.
費用はだいたい 10,000-15,000 円/月くらいのようです.
東京以外にもあるかもしれませんが,
ちょっとそこまでは調べきれていません.
それ以外で知ってるところというと,
「数学カフェ」などの勉強会をやっている人達がいます.
そうしたところでやるのも一手でしょうね.
例えば次のところはちょっとやりとりしたことがあります.
東京と埼玉でやっていますね.
これは会場代やお茶代くらいで「有料」という感じではないようです.
その他には Skype を使ってオンラインセミナーをしている人達がいたりもします.
私も以前, 東大と京大の学生がメインの参加者だった Skype セミナーに参加したことがあります.
オンラインならいろいろやりようもあるし,
私にも多少のノウハウがあるのでそれをご案内することもできます.
目的はいろいろあるでしょう.
昔の苦手意識を克服したいだとか,
大学の数学のようにもっと進んだことにも興味はあるけれど,
まずは中高数学をちゃんとやりたいとか.
一応書いておくと,
中高数学を知っているかどうかと大学の数学への適性・耐性があるかは全く別の問題です.
さらに物理なり化学なり生物なりの理工学をやるのに,
中高の数学を知っていることがどの程度意味があるかも微妙なところです.
必要な数学は大学教養の数学みたいなところでもっとハードですし,
1 から勉強し直しという感じになるからですね.
全く無駄とまではいいませんが.
で, 中高数学の復習がしたい,
そういうニーズに応えるには圧倒的なマンパワーがいります.
「どこがわからないのかがわからない」みたいなことを言う人もいるからですね.
たいてい何 1 つわかっていない状態です.
これをきちんと納得してもらった上で必要なところに戻らなければいけなくて,
ものすごい時間とパワーが取られます.
あなた 1 人で対処できるならそもそもそんな状態になっていないはずで.
指導者がちゃんとマンツーマンで付き合って,
もつれた糸を解きほぐす必要があります.
そういう意味でもたくさん指導者がいて選ぶことができる,
塾のようなふつうの中高生向けサービスを使うのがお勧めです.
地域ごとにも特色があるようで何がいいかは個別にきちんと調べる必要があります.
大手だとどうしても画一的な対応になりがちなので,
個人経営レベルの塾がよさそうです.
大人相手でもじっくり付き合ってくれるのは,
個別に小回りが効くところに行った方が早いんじゃないかと.
もちろんあなたの身近にリアルなサービスがないかもしれません.
ペースメーカーという意味では既存の通信教育は使えるでしょう.
しかしその性質上, 中学高校の勉強の文脈,
もっとはっきり言えば受験対策が基本です.
もちろんそれで良ければ問題ありません.
それで駄目なら, 例えばその数学がどんな役に立つか知りたいと思うなら,
その手のサービスだと明らかに不十分ですね.
「こんなことが知りたい」と私に連絡してもらえれば,
私の知る限りで本などの適当なコンテンツを紹介できます.
でもこの本を 1 人で読み切れない問題が解消できません.
ペースメイクも厳しいですね.
それにその手のコンテンツは「帯に短し襷に長し」で,
指導者なしでの扱いがけっこう難しいです.
だからこそこの講座を作ったわけですし.
必要ならお問い合わせフォームなどから連絡してください.
ここまで書いた上であえて勧めるなら,
やはり先程紹介した数学ガールですね.
シリーズ揃えるとそれなりのボリュームになるのは難点と言えば難点です.
ただこのシリーズがいいのは,
登場人物が一緒に悩んでくれることです.
よく対話形式の参考書もありますが,
結局生徒サイドも相当要領がよくて「そんなにすぐわかるか」と言いたくなることが多いです.
その点, 数学ガールの登場人物は割と物分かりが悪いし,
突っ込みどころにガンガン突っ込んでいくので,
ふつうの中高数学の復習や再入門とは一味違った楽しさがあります.
この節はある程度中高数学をちゃんとできている前提で書きます.
そうでないと書きづらいですし.
何はともあれまずは 2 冊本を紹介しておきます.
両方とも今回のメインテーマ, 微分方程式に関する本です.
それぞれ詳しい書評を次のページに書いています.
書評ページでは本の記述にさらに踏み込んだことまで記録しています.
ぜひ参考にしてください.
今回の講座で設定した数学的なレベルから見て接続がいいのは,
後者の『微分方程式の定式化と解法』です.
数学的に踏み込んだ面白さとしては前者の『常微分方程式の新しい教科書』です.
具体的なプログラミングに使うのは難しいですが,
シミュレーションを含めた数値計算に関しては次の 2 冊を勧めます.
前者は中学高校でやってきたことをどうやって計算機に計算させるか,
それを詳しく議論しています.
もちろん微分方程式の話もありますよ.
後者は計算機で計算する, つまりプログラミングして計算するときにどんな注意が必要かを解説しています.
もしあなたがコンピュータの計算は厳密なんだと思っているなら,
衝撃を受けるかもしれません.
私は高校生のとき, 東大の計数工学科 (当時) のオープンキャンパスに行ったときにこうした話を聞いて,
驚いたことがあります.
この講座の続きとしてはまずこの辺の本をお勧めしておきます.
歯ごたえのあるラインナップです.
中高生向けの微分方程式の本というのもなかなかないので,
厳しいところですね.
上で紹介した本や微分方程式以外の方向性についてもう少し.
大雑把に言って数学科の数学方面か,
物理なり何なりの「応用」方面かを想定しています.
あなたがやりたいかどうかにもよりますが,
理工系の基礎としてやはり物理があるので物理は割と誰もが通る道です.
物理向けの数学という視点からはいろぶつこと琉球大学の前野昌弘さんによる次の本があります.
ヴィジュアルガイドの名の通り,
図がたくさん使われています.
前野さんじたい物理学者なのでその観点から見た解説です.
実際に物理では多変数を扱う必要があるので,
きちんとやるなら明らかに不足はありますが,
中高の数学からの接続という点ではむしろいいところでしょう.
あなたは文系かもしれないので文系向け数学,
特に統計学という方面もありますね.
残念ながらというか,
文系でも数学を使わなければいけない分野があります.
経済でも微分方程式が出てくること,
本編で紹介しましたよね?
さらに言うと何かのデータを処理しないといけないなら,
その時点で統計学が出てきます.
そこでは微分積分や指数関数の処理が必要なので諦めてくださいね,
と言わざるをえません.
最近だと文学でも文体研究でテキストをデータにして,
自然言語処理なりの統計処理をかけた結果を使って研究する話もあるようですし,
文学部でも使う人は使うでしょう.
自然言語処理は携帯の漢字かな変換のようなところで使う技術ですね.
プログラミングも必要です.
個別具体的な本もいろいろ知ってはいますが,
一般的には書きづらいところがあります.
実際にもっと専門的な内容の無料の通信講座,
現代数学観光ツアーに参加された方で,
ガロア理論のような大学数学に興味があるものの,
中高の数学もままならないという方がいました.
中高の数学を知っているからと言ってガロア理論が勉強できるわけでもありません.
ただ現状を詳しく把握できないとどこからどう勧めたらいいのか,
何とも言えないところがあります.
人に合わせて興味があるところからやっていくのがやはりよくて,
そこをきちんと擦り合わせないと勉強が
不当につらくなってしまいます.
あなたはプログラミングに興味があるかもしれません.
言語から何からいろいろな観点があるのですが,
グラフを描こう, 科学技術計算をしようという観点からは Python がいいと思っています.
実際この講座で紹介したプログラミングのコードは Python のコードです.
今勉強を兼ねてプログラミング関係のコンテンツや記事も書いているところです.
まとまった形にはしていないのですが,
例えば私のサイトでプログラミングや Python のカテゴリを見てもらえれば,
記事がたくさん置いてあります.
適宜参考にしてください.
この記事の中にも Python 入門といった記事はあります.
ただ他のサイトのコンテンツの方が網羅的で取り組みやすいのが正直なところです.
まだそこまできちんと整備していないので.
少し前まではドットインストールをお勧めしていたのですが,
全部勉強するのに有料会員になる必要が出てきてしまいました.
こちらも少しずつ情報を整備していく予定です.
漫画や小説からノンフィクション,
数学者のエッセイまでいろいろ入れてあります.
数学的に異常な難易度を誇る本も入れてあるので注意してください.
実際に私が息抜きに気軽に適当に眺めている本で,
基本的にどれもお勧めできるクオリティです.
今回もアンケートがあります.
改善につなげるためぜひ回答をお願いします.
この講座はいったんここで終了です.
ここまでお疲れ様でした.
細かい部分をあまり説明していませんし,
けっこう大変な内容だったと思います.
今後も数学や物理の学習に関する情報は配信していきます.
ぜひお付き合いください.
\clearpage
こちらに PDF があります.
サイトでは見づらい方, コンテンツを手元に置いておきたい方は
ダウンロードしてご覧ください.
前回数列の説明をしました.
関数は数と数の対応ルールのことで,
数列は関数の特別な形,
つまり自然数と実数の対応ルールのことでした.
で, 漸化式.
漸化式は数列のルールの作り方の 1 つなのでした.
ここから微分方程式から見た微分の発想につながっていきます.
この添字, 普段は数列の $n$ 番目としか言ってないですね.
もっと積極的な意味をつけましょう.
それは$n$ 分後とか$n$ ステップ後とか時間的な意味です.
2 つめの $n$ ステップがわかりづらいでしょう.
それはあとで説明します.
小学校の頃こんな問題があったのを覚えているでしょうか?
A さんは毎分 75 m の速さで歩いて家を出ました.
A さんは 5 分後に何 m 先にいるでしょうか?
こんなのは $5 \times 75 = 375$ で一発です.
これを数列のスタンスで解いてみましょう.
A さんの $n$ 分後の家からの距離を $a_{n}$ と書くことにします.
まず $a_{0} = 0$ ですね.
毎分 75 m の速さで進むので 1 分後の家からの距離 $a_{1}$ は $a_{1} = a_{0} + 75 = 75$ です.
2 分後の距離 $a_{2}$ は 1 分後の距離から 75 m 追加です.
式で書けばもちろん $a_{2} = a_{1} + 75 = 150$ ですね.
これをくり返せば $n+1$ 分後の家から距離 $a_{n+1}$ は $a_{n+1} = a_{n} + 75$ です.
1 分前にどこにいたかわかれば次にどこにいるかわかります.
いま毎分 75 m と言ったから $n$ 分後としただけで,
これが毎時 75 m だったら $n$ 時間後を考えたいし,
毎秒 75 m だったら $n$ 秒後を考えたいですね?
いちいち状況に合わせて $n$ の呼び名を変えるの面倒ですから,
$n$ は時間というよりこの適当な単位のことだと思いたいです.
それを言葉ではっきりさせるために $a_{n}$ を $n$ ステップ目と呼びましょう.
あなたは当たり前のことを言っているだけだと思ったかもしれません.
その通りです.
そしてこの当たり前を過激に推し進めると微分に近づいていきます.
高校の力学の基本は等加速度運動です.
つまりある時刻を基準にしてその $n$ ステップ後の速さを $v_{n}$ としましょう.
速度が 1 ステップごとに $\alpha$ だけ増えるのなら,
等加速度運動は $v_{n+1} = v_{n} + \alpha$ で決まる漸化式で書けます.
実際には 1 ステップとかケチくさいこと言わないで一気に次のように書きます.
基準時刻から任意の時刻 $t$ が経過したあとの速さは $v(t) = v_0 + \alpha t$.
ここで $v_0$ は初期速度である.
漸化式を解くと $v_{n} = v_{0} + \alpha n$ なので,
形式的に $n$ を $t$ に変えればいいだけです.
でもいきなり上みたいに書かれて物理全然わからん!と思ってたりしませんでした?
やってること, 小学校の頃と何も変わりませんよ?
ふつうは適当に平均を取ると毎分 75 mなのであって,
ずっと同じペースで動いているわけないですね.
漸化式で言うと $v_{n+1} = v_{n} + \alpha$ みたいなのが成り立ちません.
$\alpha$ がステップごとに変わるわけで.
どこにどう責任を押しつけるかは難しいところです.
いちばんシンプルなのは上に書いたように $\alpha$ を $n$ の関数だと思うことですね.
つまり $v_{n+1} = v_{n} + \alpha_{n}$ と $\alpha$ を定数から関数 (数列) にしてしまいます.
面倒なのでもう考える問題は A さんが場所 B から場所 C に移動する状況にしましょう.
そして $\alpha_{n}$ の部分をなるべくシンプルにしたいです.
実はここで究極的な理想形を考えると微分が出てきます.
とてもシンプルです.
ただもう少し現実の泥くさいところ・面倒くさいところを考えていきます.
何をどうがんばっても $\alpha_n = \alpha$ にすることだけはできません.
逆にもしそうできてしまったとすると全ての運動が等加速度運動にしかならないからです.
残念ながらそんな世界には生きていません.
とにかく $\alpha_{n}$ が曲者です.
これを何とかしたい.
もちろん一般的にはどうにもなりません.
でも特殊な状況なら何とかなるかもしれません.
数列と漸化式を基本に考えていても苦しいので,
いったんそこから抜けることにしましょう.
最初に漸化式は微分方程式を近似して出てくることを見ましたね?
このミニ講座では微分方程式が本体です.
何の指針もなく五里霧中を彷徨うよりも,
大事な微分方程式と関係した漸化式を考える方が手が打ちやすいです.
というわけで数列のレベルでぐちゃぐちゃ考えるのはこの辺にして,
もう微分に行ってしまいましょう.
今回はこの辺で終わりにします.
お疲れ様でした.
今回もアンケートがあります.
改善に役立てたいのでぜひご協力をお願いします.
ではまた次回をお楽しみに!
リストの一番左は特別: そこには手続きを置く.
命名と抽象化について割といいことが書いてある気がする.
太字で強調してあるところがあって,
何となく大事そうだからこちらにも書かせてもらおう.
プログラミングは命令的知識を扱い, 数学は宣言的知識を扱う.
命令的というのは「如何にして為すか」という意味で,
宣言的というのは「何であるか」という意味.
素数夜曲では「函数」と書いているがめんどいので関数と書くことにする.
いわゆるラムダの話から入る.
次の関数は $f(x) = x^2$ だ.
ラムダで書くと $\lambda x.x^2$.
ラムダ以外の書き方もあるがとりあえず本の順に沿って書いていこう.
(lambda (x) (* x x))
(write ((lambda (x) (* x x)) 5))
RESULT
25
一次関数 $x \mapsto ax + b$ は Scheme だと次のように書く.
(lambda (a) (lambda (b) (lambda (x) (+ (* a x) b))))
(write ((((lambda (a) (lambda (b) (lambda (x)
(+ (* a x) b))))
2) 3) 5))
RESULT
13
こんな鬱陶しいのは嫌だ.
LISP が嫌いになる理由として括弧が挙げられる理由を強く感じる.
慣れるとむしろ便利なくらいと聞くが道は遠い.
何はともあれもう少し簡単な書き方がラムダレベルでもあるし Scheme にもある.
$\lambda abx.ax+b$ で, Scheme は次の通り.
((lambda (a b x)
(+ (* a x) b))
2 3 5)
(write ((lambda (a b x)
(+ (* a x) b))
2 3 5))
RESULT
13
素数夜曲でいつ出てくるのかわからないが,
大分長いこと lambda をひっぱるようだ.
鬱陶しいので適当に検索してきて lambda 抜きの関数定義を書いておく.
(define (hello name)
(string-append "Hello " name "!"))
(write (hello "World"))
RESULT
"Hello World!"
$f(x) = x^2$ を書いてみよう.
myfunc = lambda x: x ** 2
print(myfunc(2))
RESULT
4
普通の関数で書いてみよう.
def f(x):
return x**2
print(f(2))
RESULT
4
Haskell での lambda は次の通り.
mysquareInt = \x -> x ^ 2
mysquareDouble = \x -> x ** 2
main = do
print $ mysquareInt 2
print $ mysquareDouble 2.0
RESULT
4
4.0
一般的な lambda と無名関数の違いがよくわかっていない.
var f = function(a) { return Math.pow(a, 2) };
console.log(f(3));
var a = (function(a){ return Math.pow(a, 2) })(3);
console.log(a);
console.log((function(a){ return Math.pow(a, 2) })(3));
RESULT
9
9
9
いわゆる束縛変数の文字を変えるというやつ.
Lambda calculus だし,
たぶんいわゆるラムダ計算なのだろう.
プログラムじたいはないのでばっさり省略.
省略.
略.
略.
略.
対象に名前をつけてプログラム上のどこからでも使えるようにすることを,
トップレベル定義 (top-level definition) という.
1 次関数を定義する.
(define linear
(lambda (a b x)
(+ (* a x) b)))
(write linear)
(newline)
(write (linear 2 3 5))
RESULT
#<procedure linear (a b x)>
13
最後の (linear 2 3 5) は $2 \times 5 + 3$ を計算している.
Lambda を使わない省略記法もあって MIT 記法 という.
(define (linear a b x)
(+ (* a x) b))
(write linear)
(newline)
(write (linear 2 3 5))
RESULT
#<procedure linear (a b x)>
13
この手の省略記法は糖衣構文 (syntax sugar) と呼ばれる.
トップレベルは対話的入力時のプロンプトに象徴される領域で,
大域環境 (global environment) と呼び,
この環境下で定義された変数を大域変数 (global variable) と呼ぶ.
これと対になるのが局所環境 (local environment),
局所変数 (local variable) だ.
アルファ変換した linear# も書いておこう.
(define (linear# u v w)
(+ (* u w) v))
(write linear#)
(newline)
(write (linear# 2 3 5))
RESULT
#<procedure #{linear#}# (u v w)>
13
アルファ変換は局所変数を別の文字 (列) に置き換えてもいいとか,
とりあえずそういう雑な理解をしている.
計算機科学的にもっと面倒なところまでカバーしている可能性があり,
そこまで調べていないからいまの理解の雑さは割と真剣に警戒している.
用語の定義があるから適当に書いておこう.
対象を挿入する場所 (スロット) は本文と同じくを角括弧 (<>) で書きたいが,
いろいろな事情から隅括弧 ([]) で書くことにする1.
Lambda 抜きの糖衣構文を使った define の一般形は次の通り.
定義方法: (define ([name] [fps]) [body])
利用方法: ([name] [aps])
本の記述が死ぬほどわかりづらいが,
とりあえず引用しよう.
手続の作用の対象となるのが引数 (parameter, argument) だ.
[fps] は関数内部で使われる仮引数 (formal parameters) の略記で,
[aps] は実際の処理対象である実引数 (actual parameters) の略記.
[body] は任意個数の式 (expression) で構成される手続の本体だ.
Scheme での式はアトムやリテラルを含めて評価値が戻ってくるもの全てを指す.
とりあえずここまでくり返し書いてきた関数定義を.
(define (linear a b x)
(+ (* a x) b))
linear が [name], [fps] は a, b, x で,
[body] が (+ (* a x) b) だ.
[name] は要は関数の名前,
[fps] は定義した関数の引数,
[body] は関数で実際にやる処理のこと.
よくプログラミングのマニュアルでは上のタイプのよくわからない説明書きがある.
こういうやつ.
array explode ( string $delimiter , string $string [, int $limit = PHP_INT_MAX ] )
こういうのを見たらとりあえず実際にコード片を書いて実行して確認してみた方がいい.
めちゃくちゃ引数がたくさんあってそんなことやっていられないこともよくあるけれども.
Scheme の一般評価規則にしたがわない対象,
つまり特定の対象の評価値を求めないものを特殊形式 (special form),
あるは構文 (syntax) と呼ぶ.
スペシャルフォームと片仮名書きされているのをよく見かける.
ここからは代表的な特殊形式を紹介していくようだ.
一応改めて Python や Haskell での関数定義を書いておこう.
まずは変数定義.
a = 1
print(a)
RESULT
1
数を 2 乗する関数を定義する.
def f(x):
return x ** 2
print(f(2))
print(f(3))
RESULT
4
9
私が理解している限り, いわゆる変数はない.
定数関数があってそれが定数の役割をしてくれるという理解.
two :: Int
two = 2
main = do
print two
RESULT
2
こちらは整数を 2 乗する関数を定義する.
mySquare :: Int -> Int
mySquare x = x ^ 2
main = do
print $ mySquare 2
print $ mySquare 3
RESULT
4
9
use strict がないバージョン.
グローバル変数になるのでよくないと言われるやつ.
a = 2;
console.log(a);
RESULT
2
use strict はどこまで使えるようになっているのだろうか.
use strict をつけると上のキーワードなしの宣言でエラーになる.
次のキーワードはスコープとかの問題で詳しくはこことかで適当に調べる.
"use strict";
var a0 = 1;
console.log(a0);
var a = 2;
console.log(a);
let b = 3;
console.log(b);
const c = 4;
console.log(c);
RESULT
1
2
3
4
proc: procedureexp: expressionnum: numberpredi: predicatefunc: functionvar: variableint: integerconsq: consequentargs: argumentval: valuestr: stringaltna: alternativechar: characterobj: objectlst: liststm: stream変数定義的な意味での define の使い方.
正確には次のように言うべきのようだ.
記憶領域のある場所に
aという名前を与え, その場所と値 3.14 を結びつける.
値を参照するラベルが a とかいうやつだろう.
(define a 3.14)
(write a)
RESULT
3.14
define の次の用法の例でもある模様.
定義方法: (define [name] [exp])
次のように読めばいいようだ.
[name] とは [exp] の名前である.
a を定義したあとなら a を入力してもエラーが出ない.
REPL 的なアレで挙動確認しているわけではないので,
次のコード片では write を使っている.
(define a 3.14)
(write a)
RESULT
3.14
このように変数に束縛された値を確認する方法を変数参照 (variable reference) という.
次の実行結果を見ると a と定数 3.14 が関連づけられていることがわかる.
(define a 3.14)
(write a)
(newline)
(write (* 2 a))
RESULT
3.14
6.28
この a に値を格納するのを束縛 (bind) と呼ぶ.
この束縛情報の全体を環境 (environment) と呼ぶ.
定義方法: (lambda ([fps]) [body])
利用方法: ((lambda ([fps]) [body]) [aps])
改めてコード例を書いておこう.
数を 2 乗する関数を定義する.
(lambda (x) (* x x))
(write ((lambda (x) (* x x)) 2))
RESULT
4
関数定義の実験として前置記法 (prefix notation) を中置記法や後置記法で書いている.
前置記法について prefix を作っている.
(define prefix
(lambda (proc a b)
(proc a b)))
(write (prefix + 2 3))
(newline)
(write (prefix * 2 3))
RESULT
5
6
これを中置記法 (infix notation), 後置記法 (postfix notation) で書く.
(define infix
(lambda (a proc b)
(proc a b)))
(write (infix 2 + 3))
(newline)
RESULT
5
(define postfix
(lambda (a b proc)
(proc a b)))
(write (postfix 2 3 +))
RESULT
5
Python, Haskell, JavaScript のコード例はいらないだろう.
手続の評価値は内部の最後の式で決まる.
(define (linear x)
(define a 2)
(define b 3)
(+ (* a x) b))
(write (linear 5))
(newline)
RESULT
13
次の式では最後の (* a b) の結果が返る.
(define (linear1 x)
(define a 2)
(define b 3)
(+ (* a x) b)
(* a b))
(write (linear1 3))
RESULT
6
関数内部で定義した変数は外で参照できない.
a は関数の中で定義されているだけでトップレベルで定義されていないから,
関数の外で呼ぼうとするとエラーになる.
次のコードを org-babel で実行するとエラーになる.
よくわからないが guile -s で実行した結果を貼っておこう.
(define a 3.14)
(write a)
(newline)
(define (linear2 x)
(define a 2)
(define b 3)
(+ (* a x) b)
(* a b))
(write (linear2 3))
(newline)
(write a)
RESULT
An error occurred.
guile -s での実行結果
3.14
6
3.14
同じようなスコープの問題がある.
JavaScript だとグローバル変数に関する有名な問題がある.
昔は必ず var を使えというのがあった.
最近は var よりも use strict と let, const を使うのがいいのだろうか.
識者のご意見求む.
lambda の持つ列挙機能を含む特殊形式が begin.
式を対象にする限りは同じらしい.
begin だとトップレベル定義になる.
(begin
(define a 2)
(define b 3))
(write a)
(newline)
(write b)
上のコード, org-babel で実行できない.
guile -s での実行結果を書いておこう.
guile -s での実行結果
2
3
Python, Haskell, JavaScript で似たようなのがあるだろうか?
今回はこの辺で終わろう.
数学・物理学習とプログラミングと絡められないかと試行錯誤している.
数学からの取っつきやすさ,
物理からの取っつきやすさ,
プログラミングからの取っつきやすさをそれぞれ考えないといけない.
ある程度の体系性もほしい.
数学や物理としても面白い内容にしたいし,
プログラミングから見ても意味がある内容にしたい.
どうしたものかとずっと思っていたのだが,
素数という数学のキラーコンテンツとプログラミングを絡めた素数夜曲があった.
微妙なところも多いとは思いつつ,
メインエディタとして Emacs を使っているし,
LISP 系の言語はきちんとやりたいとも思っている.
そんなところでいい感じの落とし所という気がしたので,
少しずつ読み進めつつコードの記録をしていこうと決めた.
数学・物理の数値計算との相性の良さ,
私の勉強してみたさとも合わせて Python や Haskell のコードも書いていけたらいいな,
と思っている.
あと何となく JavaScript もやってみよう.
JavaScript は動きが速すぎて何かやってもすぐに動かなくなりそうで嫌なのだが,
ここでやるくらいのことならそう簡単に陳腐化しないだろうと勝手に信じてやっていこう.
基本的には付録のプログラミングパート,
特に付録 B から記録をつけていく予定だ.
Emacs の org-babel で書いているので,
素数夜曲本編とは見た目ちょっと違うコードを書いていくかもしれない.
Scheme, org-babel ともにあまりよくわかっていないがとりあえず進める.
org-babel の Scheme が標準で Guile を使っていて,
その変更の仕方がわからなかったのでとりあえず Guile を使う.
これまで入れていた Gauche を使いたかったが org-babel が対応していないようだ.
そうでないなら MIT/GNU Scheme が良かった気もするがこれも設定がわからない.
でははじめよう.
(write "write test")
(newline)
(display "display test")
RESULT
"write test"
display test
とりあえず Hello, World! で.
(write "Hello, Scheme!")
RESULT
"Hello, Scheme!"
加減乗除は次の通り.
前置記法なのがポイント.
割り算で分数にしてくれるのは割とポイント高い気がする.
(write (+ 5 3))
(newline)
(write (- 5 3))
(newline)
(write (* 5 3))
(newline)
(write (/ 5 3))
RESULT
8
2
15
5/3
print(5 + 3)
print(-5 + 3)
print(5 * 3)
print(5 / 3)
RESULT
8
-2
15
1.6666666666666667
まだ (?) org-babel との相性がよくないっぽい.
ここ を参考にちょっと改造.
Haskell でのあまりの計算には mod を使う.
main = do
print $ 5 + 3
print $ -5 + 3
print $ 5 + 3
print $ 5 / 3
print $ mod 5 3
print $ 5 `mod` 3
RESULT
8
-2
8
1.6666666666666667
2
2
console.log(5 + 3);
console.log((-5 + 3));
console.log(5 * 3);
console.log(5 / 3);
RESULT
8
-2
15
1.6666666666666667
例えば次の式.
\begin{align}
2 \div 3 + 5 \times 7 – 11
\end{align}
Scheme (LISP 系言語) だと括弧で細かく区切っていくから,
演算の優先度みたいな面倒なことをあまり考えなくても済む.
(define a (- (+ (/ 2 3) (* 5 7)) 11))
(write a)
RESULT
74/3
Python だと普通に書く.
a = 2 / 3 + 5 * 7 - 11
print(a)
RESULT
24.666666666666664
Python っぽくも Scheme っぽくも書ける.
演算子を括弧でくくると Scheme のように演算子を前に置けるから.
main = do
print $ 2 / 3 + 5 * 7 - 11
print $ (-) ((+) ((/) 2 3) ((*) 5 7)) 11
RESULT
24.666666666666664
24.666666666666664
console.log(2 / 3 + 5 * 7 - 11);
RESULT
24.666666666666664
素数夜曲 P.389 によると Scheme で sqrt は単項演算子らしい.
(write (sqrt 2))
RESULT
1.4142135623730951
入れ子にすれば複数回適用できる.
Scheme は関数合成あるのだろうか.
ここを見ると Gauche ならあるようだが.
(write (sqrt (sqrt 4)))
RESULT
1.4142135623730951
Python 版は次の通り.
Python は math をインポートしないといけないのがめんどい.
Python 3.5.1 を使っているが関数合成はないのだろうか.
ちょっと調べたところでは標準ではなさそうだった.
import math
print(math.sqrt(2))
print(math.sqrt(math.sqrt(4)))
RESULT
1.4142135623730951
1.4142135623730951
Haskell は関数合成 (.) がある.
main = do
print $ sqrt 2
print $ sqrt $ sqrt 4
print $ (sqrt . sqrt) 4
RESULT
1.4142135623730951
1.4142135623730951
1.4142135623730951
console.log(Math.sqrt(2));
console.log(Math.sqrt(Math.sqrt(4)));
RESULT
1.4142135623730951
1.4142135623730951
Scheme では exp を使えば出せる.
$e$ 自体は定義されていない?
(write (exp 1))
RESULT
2.718281828459045
Python の場合はやはり math モジュールにある.
import math
print(math.e)
RESULT
2.718281828459045
NumPy にもある.
import numpy as np
print(np.e)
RESULT
2.718281828459045
$e$ 自体は定義されていないので Scheme と同じく指数関数から近似値を出す.
main = do
print $ exp 1
RESULT
2.718281828459045
Math.exp() が自然対数の底による指数関数で,
Math.pow() は一般の底と指数の関数.
console.log(Math.exp(1));
console.log(Math.pow(2, 3));
RESULT
2.718281828459045
8
まずは等号 = を.
(write (= 1 1))
(newline)
(write (= 1 2))
RESULT
#t
#f
ここで #t は true, #f は false を意味している.
== ではないのかとちょっと驚いた.
ちなみに引数がたくさんある場合にも適用できる.
(write (= 1 2 3 4 5 6))
(newline)
(write (= 1 1 1))
(newline)
(write (= 1 1 1 2))
RESULT
#f
#t
#f
不等号 > >= などもある.
不等号はそれぞれ隣の項に対して適用させていった結果を出力するようだ.
(write (< 1 2 2 3))
(newline)
(write (< 1 2 3 4))
RESULT
#f
#t
等しくないことの判定は != だ.
print(1 == 2)
print(1 == 1)
print(1 != 1)
print(1 < 2)
RESULT
False
True
False
True
こちらは = だと代入になってしまうので == と重ねる.
Scheme の (= 1 2 3 4) みたいなのはどう書けばいいのだろう.
foldr で書けると思ったので foldr (==) 1 [1, 2, 3, 4] と書いてみたら怒られた.
相談したら「それ型エラー起こしてるから. ちゃんとして」と言われた.
Haskell, ちょっとしたところですぐこけるのでとてもつらい.
結論からいうと適当に 2 つリストを作り,
その要素を比較していく形で対処するのが普通らしい.
他のところでもそういう処理をするのがいいというコメントを頂いたことがあるので,
道具箱におさめておこう.
main = do
print $ 1 == 1
print $ 1 /= 1
print $ 1 < 2
let xs = [1, 2, 3, 4]
print $ and $ zipWith (==) xs (tail xs)
let ys = [1, 2, 3, 4]
print $ and $ zipWith (<) ys (tail ys)
-- print $ foldr (==) 1 [1, 2, 3, 4] -- エラーになった
RESULT
True
False
True
False
True
let を使うのが癪なので次回やる予定の lambda を使って書いてみる.
main = do
print $ (\xs -> and $ zipWith (==) xs (tail xs)) [1, 2, 3, 4]
print $ (\xs -> and $ zipWith (<) xs (tail xs)) [1, 2, 3, 4]
RESULT
False
True
見づらい気もするので無理せず関数定義した方がいい気もする.
mycompare :: [Int] -> Bool
mycompare xs = and $ zipWith (<) xs (tail xs)
main = do
print $ mycompare [1, 2, 3, 4]
RESULT
True
またコメントを頂いてしまった.
これ.
and $ ap (zipWith (<)) tail [1..4]
True
and $ ap (zipWith (<)) tail [1,2,2,3]
False
実行してみよう.
しばらくはまったが何となく ap を使うためにインポートが必要っぽい.
import Control.Monad
main = do
print $ and $ ap (zipWith (<)) tail [1..4]
print $ and $ ap (zipWith (<)) tail [1,2,2,3]
RESULT
True
False
最高か.
そして「素人はこんなこともわからない」というのを思い出せさてくれるし,
とてもいい体験だった.
ちょっと調べはしたが ap が何なのかあまりよくわかっていない.
とりあえずこれ.
myap f g = \x -> f x (g x)
main = do
print $ and $ myap (zipWith (<)) tail [1..4]
print $ and $ myap (zipWith (<)) tail [1,2,2,3]
RESULT
True
False
stackoverflow もきちんと全部読んだわけではないし,
この雑な理解でいいとも思わないがとりあえずこれで掴まえられるところもある.
「適当にしか理解していない (そして半端な理解だと割とどうでもいいところで無駄にはまる)」ことさえ理解しておけば,
あとでいくらでも調節は効くと数学でいつも経験することをいつも通り信じて追記を終える.
比較には === と == がある.
基本的に === でやるのが安全.
console.log(1 === 2);
console.log(1 === 1);
console.log(1 !== 1);
console.log(1 < 2);
RESULT
false
true
false
true
JavaScript で複数の比較はどう書くといいだろうか.
reduce を使おうと思ったが,
それは Haskell の最初の foldr でエラーになるのと同じ理由で駄目だ.
ここを見ると zip がないようだ.
ふつうにループを回すのでは面白くない.
~~素数夜曲 P.390 の記述によると,~~
~~Scheme は原則として大文字小文字を区別しないらしいが実装により異なるとのこと.~~
いわゆる処理系で変わるとかいうやつか?
安全性を考えるなら小文字に統一するのがいいらしい.
区別しないことを前提にキャメルケースを使うのも一手とか何とか.
今回はこのくらいにしよう.
下にあるようにコメントを頂いた.
Scheme の仕様は改定を重ねられているので、どの版であるかによって挙動が異なる部分があります。 Scheme の仕様はその正式名称を略して RnRS と呼ばれていて、 n の部分に何回目の改定であるかの番号が入ります。 最新は R7RS です。
「大文字小文字を区別しない」というのは R5RS までの仕様です。 (ただ、 R5RS 処理系を名乗る処理系でも区別する処理系や、切り替えられる処理系はあります。) R6RS 以降では区別するのがデフォルトです。
有理数や無限長整数などは R6RS 以外ではオプショナルな仕様です。 (R6RS で必須になった後、 R7RS では再びオプショナルになりました。) 処理系によっては桁あふれに配慮が必要な場合もあります。 有理数や無限長整数は Guile や Gauche では問題なく使えますので気にする必要はありませんが、その他に色々なところで仕様では未規定の部分があるので処理系を乗り換える機会にはつまらないことで躓くかもしれません。
思っていたよりも面倒そうな話だった.
ただ「ハマる可能性があるかもしれない」
「ハマるところはこの辺」というのが頭の片隅にあるかないかだけでも
全く違うことはこれまでの経験でよくわかっている.
こんなありがたいことはない.
こういうコメントを頂けるのが記事にしておく醍醐味だ.
地道に続けよう.
こちらに PDF があります.
サイトでは見づらい方, コンテンツを手元に置いておきたい方は
ダウンロードしてご覧ください.
ベクトル, 関数と来て今回は数列です.
簡単に復習すると,
関数は数と数の割り当てルールのことで,
数列も関数なのでした.
いろんな数列があります.
高校でもいろいろやります.
今回の話で 1 番大事な等比数列です.
長々とやっても仕方ないので軽く復習しときましょう.
$(x_n)$ を数列とします.
高校だと {$x_n$} と中括弧で書くと思いますが大学の数学的な都合で小括弧で書くことにします.
等比数列がどんな数列だったか,
つまりどんなルールで与えられているのかというとこんなルールです:
\begin{align}
x_{n+1}
=
\alpha x_{n}.
\end{align}
$x_{n}$ が 0 のときでも成り立つ形で書きました.
$x_{n}$ が 0 ではないならもちろん $x_{n+1} / x_{n} = \alpha$ と書けます.
どんな $n$ に対しても比が一定な数列だから等比数列というのでした.
このルールを使って一般項を出します.
$x_{n+1} = \alpha x_{n}$ で,
$x_{n} = \alpha x_{n-1}$ で,
$x_{n-1} = \alpha x_{n-2}$ で,
と延々続きます.
これが終わるのは $x_1 = \alpha x_{0}$ です.
順々に代入していくと $x_{n} = x_{0} \alpha^{n}$ になりますね.
この計算, はまる人ははまります.
実際私も高校の頃, 時々わけがわからなくなることがありました.
試験で慌てていてパニクったことを思い出します.
それはそれとして,
ここで細々とした計算を詳しく解説してると本筋を見失うので省略します.
数列をどんなところでどう使うのか,
それをちゃんと見てからきちんとやった方がよさそうですしね.
これがどこで出てきたかというと微分方程式を近似した方程式で出てきたのでした.
放射性物質の崩壊で出てきたのはこんなやつです.
\begin{align}
\frac{x_{n+1} – x_{n}}{h}
=
– x_{n}.
\end{align}
これを整理すると $x_{n+1} = (1 – h) x_{n}$ で,
$1-h = \alpha$ とすればまさにさっき書いた等比数列ですね.
さっき書いたように等比数列は指数関数で書けます.
で, 元の微分方程式の解も実は指数関数で書けます.
ちょっと先走ったことも書きました.
何はともあれ, 数列は関数で, 関数は数と数の割り当てルールのことで,
その割り当てルールとして等比数列っていう割り当て方があるんだ,
ということです.
微分方程式もその割り当てルールと見なせるっていうのが今書いたことです.
微分方程式と揃えた言い方として上の式は差分方程式と言うこともあります.
ちなみに上で紹介した割り当てルールはもっと一般的にできます.
これ, 要は次の項をその前の項で決まると言っているだけなので,
$a_{n+1} = a_{n}^2$ みたいにしてもいいし,
$a_{n+1} = a_{n} + a_{n-1}$ みたいに 2 つ前と関係あるようにしてもいいです.
この一般的なやつがいわゆる漸化式です.
もちろんこれ以外にもいろいろな漸化式があります.
あとで使うのでもう 1 つ漸化式を紹介します.
そのために 2 つの数列 $(x_{n})$, $(v_{n})$ を用意します.
で, 次のような漸化式を考えます.
\begin{align}
x_{n+1} – x_{n}
=
v_{n}, \quad
v_{n+1} – v_{n}
=
– \alpha^2 x_{n}.
\end{align}
2 つの数列が絡まりあう漸化式です.
高校でやる言葉を使うなら,
互いに階差数列を取ってると思ってもいいかもしれません.
これがどこから出てくるかというと次の微分方程式からです.
\begin{align}
\frac{d^2 x(t)}{dt^2}
=
– \alpha^{2} x(t).
\end{align}
この微分方程式も物理でよく出てくる方程式で,
いわゆる単振動の運動方程式です.
さっきの漸化式を出すには $v(t) = x'(t)$ としてもとの方程式も $v’ = – \alpha^2 x$ とすればいいです.
これを差分近似すると次の式が出てきます.
\begin{align}
x_{n+1}
=
x_{n} + h v_{n}, \quad
v_{n+1}
=
v_{n} – h \alpha^{2} x_{n}.
\end{align}
$h$ が入っていて係数が微妙に違いますが,
あまり気にしないでください.
この辺の計算プログラムも準備してあります.
この記事後半のプログラムパートを確認してください.
今回もアンケートがあります.
改善につなげるためぜひ回答をお願いします.
ではまた次回をお楽しみに!
復習でもう 1 回.
一番単純でしかも実際に使われる微分方程式としてまずは 1 階の線型常微分方程式を考えよう.
ちょっと不吉な例であるが放射性物質の崩壊の方程式を紹介する.
導出をしたければちゃんと物理を勉強してもらう必要がある.
ここでは物理は省略して数学に集中する.
\begin{align}
\frac{dx}{dt} = – c x.
\end{align}
厳密解は $x = C_0 e^{-ct}$ だ.
初期値を設定すれば $C_0$ はそこから決まる.
微分を単純に離散化すると次のようになる.
\begin{align}
\frac{x_{n+1} – x_{n}}{\Delta t}
=
-c x_{n}.
\end{align}
$\Delta t$ は $h$ と書くこともある.
整理すると次の通り.
\begin{align}
x_{n+1}
=
x_{n} – c (\Delta t) x_{n}.
\end{align}
これに沿って計算したのがいわゆるオイラー法.
次のセルではこれをコードに落としている.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
def radioactive_euler(nt, init = 10):
dt = 2 / (nt - 1)
# 初期条件設定
x = np.zeros(nt)
x[0] = init
for i in range(1, nt):
x[i] = x[i-1] - c * dt * x[i-1]
# ベクトル計算で書き直したい
return x
# 近似解
c = 1
nt = 101
init = 5
x1 = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x1)
# 厳密解
t = np.linspace(0, 2, nt)
x2 = init * np.exp(- c * t)
plt.plot(t, x2)
# 凡例
plt.legend(['approximation', 'exact'])
plt.show()
# 描画
#plt.savefig(img_file)
#b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
#img_str = "
" % (b64)
#print(img_str)
RESULT
こちらは単振動.
次に 2 階の常微分方程式を紹介しよう.
高校の物理で出てくるばねの振動(単振動)がまさにこの例だ.
項を増やすと減衰振動になったり、外力をつけたりといろいろなケースがある.
まずは一番単純な式を考えよう.
\begin{align}
\frac{d^2 x}{dt^2}
=
– \omega^2 x.
\end{align}
オイラー法やルンゲ-クッタ法は
1 階の方程式に対する計算法なので直接は使えない.
今は中間処理として $v = dx/dt$ を置いて計算すればいい.
これは単なる数値計算の便法ではない.
速度の意味もあるから,
という表面的な理由ではなくもっと深く解析力学の文脈で物理としても大事な視点だ.
もっといえばシンプレクティック計算法などもっといい計算法にも発展する.
とりあえずオイラー法で計算したい.
まずは微分方程式自体を書き直す.
\begin{align}
\frac{dx}{dt}
=
v, \quad
\frac{dv}{dt}
=
– \omega^2 x.
\end{align}
これをオイラー法で近似しよう.
\begin{align}
x_{n+1}
=
x_{n} + h v_{n}, \quad
v_{n+1}
=
v_{n} – h \omega^{2} x_{n}.
\end{align}
オイラー法をコードに落とす.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
def harmonic_euler(nt, init = (5, 0)):
dt = t_range / (nt - 1)
# 初期条件設定
x = np.zeros(nt)
v = np.zeros(nt)
x[0] = init[0]
v[0] = init[1]
for i in range(1, nt):
x[i] = x[i-1] + dt * v[i-1]
v[i] = v[i-1] - dt * (omega ** 2) * x[i-1]
# ベクトル計算で書き直したい
return (x, v)
omega = 2 * np.pi
nt = 101
t_range = 2
init = (5, 0)
harm = harmonic_euler(nt, init)
t = np.linspace(0, 2, nt)
# 厳密解
x_exact = init[0] * np.cos(- omega * t)
v_exact = - omega * init[0] * np.sin(omega * t)
# グラフ描画
plt.subplot(3, 1, 1)
plt.title('x-t graph')
plt.plot(np.linspace(0, 2, nt), harm[0])
plt.plot(t, x_exact)
plt.legend(['x approximation', 'x exact'])
plt.subplot(3, 1, 2)
plt.title('v-t graph')
plt.plot(np.linspace(0, 2, nt), harm[1])
plt.plot(t, v_exact)
plt.legend(['v approximation', 'v exact'])
plt.subplot(3, 1, 3)
plt.title('phase space')
plt.plot(harm[0], harm[1])
# 描画
plt.tight_layout()
plt.show()
#plt.savefig(img_file)
#b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
#img_str = "
" % (b64)
#print(img_str)
RESULT
見ての通り時間が進むごとに誤差が大きくなる.
nt を大きくすると少しはましになる.
実際に上のコードで nt を大きくして再計算してみてほしい.
あとまずいのは phase space の図だ.
この系はエネルギーが保存する系だから相空間内の軌道が閉じてほしいのにそうなっていない.
シンプレクティックにやれば解消できるようだが, とにかくここではよろしくない.
この方程式でオイラー法はよろしくないことがわかった.
とりあえずルンゲ-クッタでやってみよう.
とりあえず近似式を書く.
\begin{align}
x_{n+1}
&=
x_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}),
\end{align}
\begin{align}
t_{n+1}
&=
t_{n} + h,
\end{align}
\begin{align}
k_{1}
&=
v_{n},
\end{align}
\begin{align}
k_{2}
&=
v_{n} + \frac{h}{2} k_{1},
\end{align}
\begin{align}
k_{3}
&=
v_{n} + \frac{h}{2} k_{2},
\end{align}
\begin{align}
k_{4}
&=
v_{n} + h k_{3}.
\end{align}
次が $v$ の式.
\begin{align}
v_{n+1}
&=
v_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}),
\end{align}
\begin{align}
t_{n+1}
&=
t_{n} + h,
\end{align}
\begin{align}
k_{1}
&=
– \omega^2 x_{n},
\end{align}
\begin{align}
k_{2}
&=
x_{n} – \frac{h}{2} \omega^2 k_{1},
\end{align}
\begin{align}
k_{3}
&=
x_{n} – \frac{h}{2} \omega^2 k_{2},
\end{align}
\begin{align}
k_{4}
&=
x_{n} – h \omega^2 k_{3}.
\end{align}
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
def harmonic_rk(nt, init = 10):
dt = t_range / (nt - 1)
# 初期条件設定
x = np.zeros(nt)
v = np.zeros(nt)
x[0] = init[0]
v[0] = init[1]
def fx(t, x, v):
return v
def fv(t, x, v):
return - (omega ** 2) * x
# ベクトル計算で書き直したい
for i in range(1, nt):
xk1 = fx(dt * (i - 1), x[i-1], v[i-1])
vk1 = fv(dt * (i - 1), x[i-1], v[i-1])
xk2 = fx(dt * (i - 1/2), x[i-1] + xk1 * dt / 2, v[i-1] + vk1 * dt / 2)
vk2 = fv(dt * (i - 1/2), x[i-1] + xk1 * dt / 2, v[i-1] + vk1 * dt / 2)
xk3 = fx(dt * (i - 1/2), x[i-1] + xk2 * dt / 2, v[i-1] + vk2 * dt / 2)
vk3 = fv(dt * (i - 1/2), x[i-1] + xk2 * dt / 2, v[i-1] + vk2 * dt / 2)
xk4 = fx(dt * (i - 1), x[i-1] + xk3 * dt, v[i-1] + vk3 * dt)
vk4 = fv(dt * (i - 1), x[i-1] + xk3 * dt, v[i-1] + vk3 * dt)
x[i] = x[i-1] + dt / 6 * (xk1 + 2 * xk2 + 2 * xk3 + xk4)
v[i] = v[i-1] + dt / 6 * (vk1 + 2 * vk2 + 2 * vk3 + vk4)
return (x, v)
omega = 2 * np.pi
nt = 101
t_range = 2
init = (5, 0)
harm = harmonic_rk(nt, init)
t = np.linspace(0, 2, nt)
x_exact = init[0] * np.cos(- omega * t)
v_exact = - omega * init[0] * np.sin(omega * t)
plt.subplot(3, 1, 1)
plt.title('x-t graph')
plt.plot(np.linspace(0, 2, nt), harm[0])
plt.plot(t, x_exact)
plt.legend(['x approximation', 'x exact'])
plt.subplot(3, 1, 2)
plt.title('v-t graph')
plt.plot(np.linspace(0, 2, nt), harm[1])
plt.plot(t, v_exact)
plt.legend(['v approximation', 'v exact'])
plt.subplot(3, 1, 3)
plt.title('phase space')
plt.plot(harm[0], harm[1])
# 描画
plt.tight_layout()
plt.show()
#plt.savefig(img_file)
#b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
#img_str = "
" % (b64)
#print(img_str)
RESULT
今度の一致具合はなかなかよさそう.
あくまで見た目の感じではあるけれども.
良くなったのは一番下の図, 相空間軌道だ: ちゃんと閉じてくれた.
こちらに PDF があります.
サイトでは見づらい方, コンテンツを手元に置いておきたい方は
ダウンロードしてご覧ください.
ここまでは微分方程式が自然現象のシミュレーションに使えることを紹介し,
実際の微分方程式をいくつか見てきました.
実際のシミュレーションでは本当に微分方程式を厳密に解いているわけじゃなくて,
近似計算をしているんだ,
そしてそれは四則演算なんだというのを大雑把に紹介してきました.
近似といってもけっこうきちんと計算できますよ,
ということを見るためにもとにかくいくつかの実例を紹介することを優先にしてきました.
我慢強いあなたも「もういい加減きちんと数学の紹介して!」
と思ってらっしゃるはずです.
というわけで数学の紹介も簡単にやっていきましょう.
微分方程式なんだからまずは微分じゃないの?
と思う方もいらっしゃるでしょう.
でも今回の話の流れでは微分の話をする前にベクトルをやった方がスムーズなので,
ベクトルからいきます.
高校でベクトルは向きと大きさがある矢印のことと言われます.
それはそれで大事な見方ですが,
ここでは始点と終点がある矢印と思ってみてください.
始点と終点をつないだたくさんのベクトルで作った折れ線を方程式の近似解とみなすのが今回の話のキモだからです.
図のように曲線を適当な間隔でわけましょう.

わけた点に $A_0, A_1, \dots, A_n$ と順に名前をつけていきます.
そして $A_0$ と $A_1$, $A_1$ と $A_2$ というふうに始点と終点を入れかえながら隣の点どうしを結んでいきます.
$A_0$ から $A_1$ に向かって進むんだ,
という気持ちを表したいのでそれを $\overrightarrow{A_0 A_1}$ と書くことにします.
記号の上の右向き矢印で点 $A_0$ から点 $A_1$ に向かっている気分もはっきり書いています.
ここで何度も始点と終点をつないだ矢印を継ぎ足しています.
この継ぎ足しという図形の操作がベクトルの足し算です.
例を 1 つだけ書いておくと $\overrightarrow{A_0 A_1} + \overrightarrow{A_1 A_2} = \overrightarrow{A_0 A_2}$ です.
高校だとベクトルに関してもっといろいろなことが出てきます.
定数倍したり内積を取ったり長さを調べたりなどなど.
もちろんベクトルを使っていろいろやるならやっておかないと困ります.
でも微分方程式から数学を眺めてみる立場ではまずここさえ乗り越えればどうにかなります.
他のことは他のことをやるときに絡めてやってみてください.
折れ線をつないでいって曲線を近似する様子をもっときちんと見てみましょう.
具体的な曲線としては円を取り,
それを等分して折れ線近似した様子をプログラムで描きました.
次のページを見てください.
最後に近似がよくなっていくアニメーションもつけています.
基本的なプログラムはこの記事の後半にも載せています。
アニメーションは Jupyter との連携でやっているので,
この記事ではうまく動かせません.
上のリンクのページから確認してください.
今回はここまでです.
お疲れ様でした.
今回もアンケートがあります.
改善につなげるためぜひ回答をお願いします.
ではまた次回をお楽しみに!
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
div_num = 9
def draw_circle_by_polygonal_line(div_num):
thetas = np.linspace(0, 2 * np.pi, div_num)
xs = np.cos(thetas)
ys = np.sin(thetas)
plt.plot(xs, ys)
plt.gca().set_aspect('equal', adjustable='box')
円周上の点を $n$ 分割し折れ線で分点をつないでいった.
一旦個別に図示している.
for i in range(4, 12):
pnum = i - 3
plt.subplot(2, 4, pnum)
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])
plt.title("div = " + str(i-1))
draw_circle_by_polygonal_line(i)
draw_circle_by_polygonal_line(1000)
plt.tight_layout()
RESULT
9 分割までを重ねて 1 つの図にしてみた.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
for div_num in range(10):
draw_circle_by_polygonal_line(div_num)
draw_circle_by_polygonal_line(1000)
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])
plt.show()
RESULT
分割を細かくしていくことで精度が上がる様子をアニメーションで見てみよう.
下のプログラムを Jupyter 上で実行すると,
図の両端にある「+」「-」ボタンでアニメーションのスピードを変えられる.
適当に変えて遊んでみよう.
注意
次のプログラムは Jupyter 上で実行することを前提にしている.
他の実行環境で動くかどうかは確認していない.
ここを見ると Jupyter のアニメーションが見られる.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from JSAnimation.IPython_display import display_animation
def circle_data(div_num):
thetas = np.linspace(0, 2 * np.pi, div_num)
xs = np.cos(thetas)
ys = np.sin(thetas)
return (xs, ys)
circles = []
for i in range(4, 31):
circles.append(circle_data(i))
#fig = plt.figure(figsize=(10, 10));
fig = plt.figure();
ax = plt.axes(xlim=(-1.2, 1.2), ylim=(-1.2, 1.2));
line, = ax.plot([],[],lw=2);
plt.gca().set_aspect('equal', adjustable='box')
def animate(data):
x = data[0]
y = data[1]
line.set_data(x, y)
return line
anim = animation.FuncAnimation(fig, animate, frames=circles, interval=100)
display_animation(anim, default_mode='reflect')
こちらに PDF が置いてあります.
サイトでは見づらい方は PDF をご覧ください.
あなたが高校生以上なら社会の授業で次のような言葉を聞いたことがあるかもしれません.
人口は制限されなければ幾何級数的に増加するが生活資源は算術級数的にしか増加しないので生活資源は必ず不足する.
これは有名なマルサスの『人口論』で議論された話です.
幾何級数はいわゆる等比数列の和で,
算術級数は等差数列の和です.
マルサスは経済学者なので経済の話です.
経済の話でこういう数学ネタが出てくるわけですね.
ちなみに大学の経済学だと経済の話で微分積分が本当に出てきます.
ときどき経済学部の入試で数学が受験科目に入っていることがあります.
実際に大学に入ってから使うからです.
それも高校の理系のレベルを越えた数学です.
統計データをいろいろいじらないといけないので統計学が必要で,
そっちでも割といろいろ数学が必要です.
話を元に戻しましょう.
マルサスの人口論に合わせて人口の増え方を考えます.
ここではもっと一般的に生物の集団だと思いましょう.
まずは生物の種類が 1 種類だとして考えていきます.
生物っぽく人の数というよりも一般に個体数と呼ぶことにして,
時刻 $t$ での個体数を $x(t)$ とします1.
個体数の増加率を出生率の死亡率の差として定義します.
微分方程式的にいうとこれは $x'(t) / x(t)$ です.
増加率が一定値 $\alpha$ に等しいとすると $x'(t) / x(t) = \alpha$ だから次の微分方程式が出てきます.
\begin{align}
\frac{d x(t)}{dt}
=
\alpha x(t).
\end{align}
この微分方程式をマルサスモデルと言います.
これを解くと $x(t) = x(0) e^{\alpha t}$ です.
これで指数関数が出てきました.
等比数列は一般項が $a_n = 2^{n}$ のように指数関数で書けている数列だから,
これを足し上げていけばまさに幾何級数になりますね.
放射性物質の崩壊とは指数の肩の符号が変わっているだけで,
それ以外は同じ式です.
放射性物質の崩壊と生物の個体数の変化を追う方程式が同じ形をしているわけです.
シミュレーションの結果は放射性物質の崩壊のときと基本的に同じです.
いちおうきちんとやっておきましょうか.
まず微分方程式を次のように近似します.
\begin{align}\label{math_refuge00005}
\frac{x_{n+1} – x_{n}}{h}
=
\alpha x_n, \quad
x_{n+1}
=
\alpha x_n + h x_n.
\end{align}
これにしたがって解いてみましょう.
プログラムと結果の図はこのページの後半を確認してください.
もう少し現実的な設定にしてみます.
個体数 $x(t)$ が大きくなると食物が足りなかったり,
衛生状態の悪化で病気にかかりやすくなったりして増加率が減るはずです.
そこで増加率 $\alpha$ を $\alpha – \beta x$ ($\alpha$, $\beta > 0$) に変えてみます.
\begin{align}
\frac{dx}{dt}
=
(\alpha – \beta x) x.
\end{align}
これにはロジスティック方程式という名前がついています.
悪化を $\beta x$ と書くことに必然性はありません.
とりあえずやってみただけです.
解き方はさておき解は次のようになります.
あなたがもしこの関数を微分できるなら,
厳密解になることを計算して確認してみてください.
\begin{align}
x(t)
=
\frac{\alpha}{\beta} \frac{1}{1 + e^{- \alpha t}}.
\end{align}
これを近似すると次のようになります.
\begin{align}
\frac{x_{n+1} – x_{n}}{h}
=
(\alpha – \beta x_{n}) x_{n}, \quad
x_{n+1}
=
x_n + h (\alpha – \beta x_{n}) x_{n}.
\end{align}
これも厳密解と比較しながら計算してみます.
プログラムや図はこの記事の後半を参照してください.
これもよく近似できていることは認めてもらえるんじゃないでしょうか?
でもうるさい人は「近似の精度はどう決めるの?」みたいなことを考えているはずです.
もちろんいたって真っ当な指摘です.
ただけっこう難しい問題がたくさんあって,
なかなかスカっと綺麗に回答できません.
あとで少し考えることにしてここではこのまま進みます.
あなたは数値計算の精度とかそれ以前の問題として,
次のように思っているかもしれません.
もちろんこれも正しい指摘です.
例えばマルサスモデルでもロジスティック方程式でも右辺が気になります.
実際には成長しないと子作りできないですからね.
その成長が必要だという情報が方程式に盛り込まれていません.
生物で考えるなら食物連鎖があるわけで食べられて個体数が減ることだってあります.
天敵が増えたらその個体の数はあおりを受けて激減するはずです.
もちろんこんなことは折り込み済みで,
例えば遅延型微分方程式,
ロトカ-ボルテラ方程式と呼ばれる方程式が対応しています.
これを調べるのも面白いんですが,
今回はこのくらいにしておきましょう.
簡単にまとめると,
今回は近似の精度が気になるあなたのために,
実際どのくらいよく本当の解を近似できているのかを調べました.
経済学や生物のようにあんまり数学との関係がなさそうな分野と数学の関係も紹介しています.
長くなってきましたが最後にちょっと大事な話.
式 \eqref{math_refuge00005} で $\alpha = 1$, $h = 1$ とすると $x_{n+1} = 2 x_n$ になります.
この漸化式は高校でも出てくるやつで, もちろん公比 $2$ の等比数列の漸化式です.
底が $2$ か $e$ かの違いはあっても指数関数で書けることは同じです.
これはたまたまではなくて数学的に意味があることです.
$x_{n+1} – x_{n}$ のように数列の隣の項の差を差分と言うことがあります.
もちろん微分との関係を強く意識した言葉です.
数列の問題, もっと言えば漸化式の問題は微分方程式とも深い関係があります2.
あまり深入りはしませんが数学や物理に興味があるなら覚えておくといろいろ楽しいことがあります.
今回もアンケートがあります.
改善に繋げたいのでぜひ積極的に回答をお願いします.
ではまた次回をお楽しみに!
基本は放射性物質の崩壊と同じだからあっさり行こう.
\begin{align}
\frac{dx}{dt}
=
\alpha u.
\end{align}
厳密解は $x(t) = C_0 e^{\alpha t}$ だ.
初期値を設定すれば $C_0$ はそこから決まる.
微分を単純に離散化すると次のようになる.
\begin{align}
\frac{x_{n+1} – x_{n}}{h}
=
\alpha x_{n}.
\end{align}
整理すると次の通り.
\begin{align}
x_{n+1}
=
x_{n} + \alpha h x_{n}.
\end{align}
オイラー法でコードに落とそう.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
def malthus_euler(nt, init = 10):
dt = 2 / (nt - 1)
# 初期条件設定
x = np.zeros(nt)
x[0] = init
for i in range(1, nt):
x[i] = x[i-1] + c * dt * x[i-1]
return x
c = 1
init = 1
nt = 101
x_approx = malthus_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
t が大きくなると厳密解の方が少し大きくなる.
指数の肩の符号が違うだけなので,
放射性物質の崩壊で形式的に時間を負の方に伸ばせば同じようにズレが出てくるはずだ.
私の今の感覚だとこのズレが今回のように大きくなるように出るのかどうかまではわからない.
何はともあれ区間の分割数 nt を大きくしてみたのが次の結果:
具体的には 101 から 1001 にした.
ライブラリが読み込まれていたり関数 malthus_euler が定義されている前提のコード片なので,
これだけで実行できるわけではない.
IPython か Jupyter なら順にセルを読み込んでいけばそのまま実行できる.
c = 1
init = 1
nt = 1001
x_approx = malthus_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
$nt = 101$ よりも近似の精度が良くなった.
この辺は単純に振る舞ってくれるようだ.
調和振動子だとオイラー法自体がうまく動かないから注意したい.
方程式は次の通り.
\begin{align}
\frac{dx}{dt}
=
(\alpha – \beta x) x.
\end{align}
厳密解は次の通り.
\begin{align}
x(t)
=
\frac{\alpha}{\beta} \frac{1}{1 + e^{- \alpha t}}.
\end{align}
近似すると次の通り.
\begin{align}
\frac{x_{n+1} – x_{n}}{h}
=
(\alpha – \beta x_{n}) x_{n}, \quad
x_{n+1}
=
x_n + h (\alpha – \beta x_{n}) x_{n}.
\end{align}
では数値的に解いてみよう.
とりあえず $\alpha = 2$, $\beta = 1$ で計算している.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
def logistics_euler(nt, init = 10):
dt = T / (nt - 1)
# 初期条件設定
x = np.zeros(nt)
x[0] = init
# ベクトル計算で一気に計算したい
for i in range(1, nt):
x[i] = x[i-1] + dt * (alpha - beta * x[i-1]) * x[i-1]
return x
alpha = 2
beta = 1
init = 1
nt = 101
T = 5 # 時間変化を見る最大値
x_approx = logistics_euler(nt, init)
plt.plot(np.linspace(0, T, nt), x_approx)
t = np.linspace(0, T, nt)
x_exact = (alpha / beta) / (1 + np.exp(- alpha * t))
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
近似解と厳密解にちょっとズレが見える.
そこで nt = 101 から nt = 1001 にしてみた.
alpha = 2
beta = 1
init = 1
nt = 1001
T = 5 # 時間変化を見る最大値
x_approx = logistics_euler(nt, init)
plt.plot(np.linspace(0, T, nt), x_approx)
t = np.linspace(0, T, nt)
x_exact = (alpha / beta) / (1 + np.exp(- alpha * t))
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
こちらに PDF が置いてあります.
サイトでは見づらい方はこちらをご覧ください.
今回アニメーションはありませんが引き続きシミュレーションの話をします.
今回のテーマは近似についてです.
シミュレーションというのはいいし,
近似というのもいいけど,
実際どのくらいよく近似できてるの?
こういう問題を考えてみましょう.
ちゃんと使えばすごく精度はいいですよ,
というのが今回の主張です.
こんな疑問に答えるのが今回のテーマです.
で, ちょっと不吉な例ではあります.
しかしというか残念ながらというか,
役に立つ例になってしまったので放射性物質の崩壊に関する微分方程式を考えてみましょう.
いちおう微分方程式から書いておきますね.
\begin{align}\label{math_refuge00002}
\frac{dx}{dt}
=
– c x.
\end{align}
気分だけ説明しておくと放射性物質が単位時間あたりに崩壊するスピードはそのときの放射性物質の量に比例する,
というのを式で書いた形です.
$c$ には物理的に大事な意味がありますが式や計算を軽くするためにここでは $c = 1$ にして話を進めます.
そして前回のようにこの微分方程式を近似した式を書いてみます.
\begin{align}\label{math_refuge00003}
\frac{x_{n+1} – x_{n}}{h}
=
– x_{n}.
\end{align}
これは時間間隔 $h$ が十分小さければ時間 $h$ だけ離れた放射性物質の量 $x_{n+1}$ と $x_{n}$ の関係が上の式で書けることを意味しています.
物理的な話はこのくらいにして,
この式にしたがって近似計算 (シミュレーション) したとき,
本当の答えと近似した解がどのくらい近いのかを考えてみます.
何でこれにしたかというともとの方程式の答えが厳密に書けるからです.
あなたが指数関数, 特に自然対数の底 $e$ をご存知なら,
式 \eqref{math_refuge00002} の答えは $x(t) = C e^{- t}$ と書けます.
指数関数の微分をご存知ならこれを \eqref{math_refuge00002} に代入してみて本当に方程式の答えになっていることを確認してみてください.
この答えと式 \eqref{math_refuge00003} にしたがって計算した答えがどのくらい一致するか,
コンピュータに計算させてみましょう.
結果はこの記事の後半,
プログラムパートに載せてあります.
あなたは「近似なんて雑なことしていいの?」なんて思っているかもしれません.
でも記事の後半にある「比較のために重ねる」節のグラフを見るとわかるように,
厳密な様子ともよく合っているのでちゃんと使えば問題ありません.
「ちゃんと使う」の「ちゃんと」が難しいと言われればそれはもちろんそうなんですが,
それは本当に難しい話です.
せっかくなので分割をどんどん大きくしていって近似の精度が上がっていく様子もグラフにしておきました.
プログラムが書けるとこういうのもささっとやれます.
遊びやすくなるのは本当にありがたいですね.
中高生の頃にこれを知っていたら数学や物理だけじゃなくてプログラミングにもはまっていたと思います.
やってみるとわかるんですが,
書いたプログラムで意図通りの結果が出るとけっこう嬉しいです.
パラメータをいじったり方程式を変えたり,
ちょろっといじって遊べる要素も増えて遊べる範囲が広がるのもいいですね.
Python プログラミングに関する資料,
数値計算・シミュレーションに関する資料は GitHub に上げてありますし,
今後も講座の進展に合わせて少しずつ増やしていくのでぜひあなたも遊んで遊んでみてください.
このサイトにもいろいろな情報がありますよ.
あと今回の中高数学上のポイントを復習しておきます.
指数関数とか自然対数の底とか,
はたまたその微分が出てきました.
放射性物質の崩壊をちゃんと調べるにはこういう数学が必要なんです.
次回はまた別の微分方程式を紹介します.
経済学や生物学で出てくる微分方程式ですよ.
これが終わったらもう少し数学的にちゃんとした話をしましょう.
まずは微分方程式の射程距離の長さを知ってほしいからです.
もっとちゃんと数学の説明してほしいというあなた,
もうちょっと我慢してください.
今回もアンケートがあるのでぜひ回答をお願いします.
ではまた次回をお楽しみに!
一番単純でしかも実際に使われる微分方程式として,
まずは 1 階の線型常微分方程式を考えよう.
ちょっと不吉な例であるが放射性物質の崩壊の方程式を紹介する.
導出をしたければちゃんと物理を勉強してもらう必要がある.
ここでは物理は省略して数学に集中する.
\begin{align}
\frac{dx}{dt}
=
– c u.
\end{align}
厳密解は $x(t) = C_0 e^{-ct}$ だ.
初期値を設定すれば $C_0$ はそこから決まる.
微分を単純に離散化すると次のようになる.
\begin{align}
\frac{x_{n+1} – x_{n}}{\Delta t}
=
-c x_{n}.
\end{align}
時間っぽく見えるよう $\Delta t$ と書いてみた.
タイプが面倒なので $h$ と書くこともある.
整理すると次の通り.
\begin{align}
x_{n+1}
=
x_{n} – c (\Delta t) x_{n}.
\end{align}
これに沿って計算したのがいわゆるオイラー法.
コードに落としていこう.
まずは厳密解 $x(t) = e^{-t}$ をグラフに描こう.
初期値は 1, 上の定数も $c = 1$ としている.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
c = 1
init = 1
nt = 100
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['exact'])
plt.show()
RESULT
今度は近次解を計算してグラフにしてみよう.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
def radioactive_euler(nt, init = 10):
dt = 2 / (nt - 1)
# 初期条件設定
x = np.zeros(nt)
x[0] = init
for i in range(1, nt):
x[i] = x[i-1] - c * dt * x[i-1]
# ベクトル計算で書き直したい
return x
c = 1
init = 1
nt = 101
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
plt.legend(['approximation'])
plt.show()
RESULT
これだけではよくわからないので重ねて描いてみる.
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
ほぼピッタリ重なって見える.
シミュレーションの精度はけっこう良さそうだ.
もちろん細かいことを言えばいろいろあるけれども,
それはもっと進んだお話なので一旦不問にしておこう.
上のグラフは区間を 100 分割して計算した.
nt = 101 と言うところで分割を指定している.
101 と言うのは分点数の指定なので実際の分割としては 100 等分なのだ.
ここで刻みが荒いといい近似にならないことを大雑把に調べておこう.
nt = 3
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
nt = 101
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
何となくそれっぽい形にはなっている.
もちろん精度は大したことない.
nt = 5
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
nt = 101
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULTS
あまり適当なことを言うのも良くないが, けっこう滑らかに見える.
精度はまだまだ.
nt = 11
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
nt = 101
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
まだ差が目で見える.
nt = 51
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
nt = 101
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
100 は最初にやってあるから一気に大きくしてみた.
nt = 1001
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
nt = 101
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
コンピュータで計算するとき,
いろいろな都合があって分割数を多くすれば単純に精度が上がっていくと言うわけではない.
桁落ちや丸め誤差のような問題がある.
それはそれで別途検討する必要があるのだ.
こちらに PDF が置いてあります.
サイトでは見づらい方は PDF をご覧ください.
今回はイントロとして説明することがいくつかあります.
私が作っている他の講座と同じように,
この中高数学駆け込み寺では細かいことを詳しくやっていくよりも,
数学を勉強するモチベーションになるように数学の大きな姿をお見せしていきます.
で, 実際に何をするのかというと,
微分方程式のシミュレーションを通じて中高数学の大枠を掴んでいきます.
もっと具体的に言うと微分方程式のシミュレーションには中高でやる数学がほとんど全て出てきます.
使われているシーンを具体的に見てもらってモチベーションにしてもらおう,
そんな講座です.
この講座では文字計算には慣れていることを前提にしています.
あなたが中高生なら微分は聞いたこともないかもしれません.
あなたが中高の数学の復習をしたいと思っている大人なら,
「自分にはまだ微分なんて早いんじゃ?」と思っているかもしれません.
メインの微分方程式が何だかわからないんじゃ読み進めるのはきついですよね?
というわけでまずは微分方程式の説明をします.
微分方程式は理工系の基礎です.
いわゆる自然法則は微分方程式で表現されることが多いのです.
これがどこで使われるかというと,
例えば洪水が起きたときの被害予測に使われます.
どんな規模で起きた洪水が市街地のどこまで進入してくるか?
こういうシミュレーションをテレビで見たことがないでしょうか.
このシミュレーションに使われているのが微分方程式です.
他にも微分方程式はありとあらゆるところで使われています.
ゲームの CG を自然に見せるためにはまさに自然法則に則って風や水の動きを表現しないといけません.
だから微分方程式がその背後に隠れています.
天気予報はいまの気象条件から未来の様子を予測する必要があります.
この予測にも微分方程式を使っています.
機械を動かしていると熱くなることがありますね?
あまりにも熱くなりすぎると機械が暴走してしまうので熱を効率よく逃がす必要があります.
そのためには熱の流れを考えてその流体の動きをきちんと制御する必要があります.
ここでは熱の流れの微分方程式を考えてそれをシミュレーションします.
車を作るとき空気抵抗を調べるためにも流体力学が必要で,
シミュレーションも使ってモノづくりにも活かしています.
他にもいちいち挙げきれないくらい身の回りに微分方程式のシミュレーションを使っているモノがあります.
微分方程式自体をきちんと調べようと思うと大学の数学が必要です.
でもシミュレーションを中心に考えれば中高の数学でかなりいろいろなことがわかります.
わかるだけじゃなくて実際に微分方程式を解いて遊んでみることだってできます.
この自分でも遊べるところまで持っていけるのが大事じゃないかと思っていて.
それが微分方程式のシミュレーションを選んだ理由の 1 つです.
さて, ここまでで微分方程式が何で大事かはわかってもらえたでしょうか?
これをやれば数学をいろいろ遊び倒せそうと思ってもらえたら嬉しいです.
ちゃんとがんばれば自分でゲームを作ったりもできますしね.
次はもう少し数学的な内容に踏み込みましょう.
何はともあれ微分に関する話をやるわけです.
そして微分は高校でやる内容です.
あなたは微分に対して嫌な思い出を持っているかもしれません.
あなたは「中学レベルからやり直したいのにそんなの無理だ!」と思っているかもしれません.
もっと簡単なところからやってほしいと思っているかもしれません.
でもこの講座では中高数学の大枠を掴んでもらう講座です.
そして今回はこの講座の大枠を掴んでもらう講座です.
どうかもうしばらく辛抱して読み進めてください.
まず何をするのかを見てもらいたいので,
次のページを開いて 1 番下までスクロールしてください.
動画の再生ボタンがありますから動画を再生してみましょう.
これは 1 次元移流方程式のシミュレーション結果です.
プログラムを書いてコンピュータに計算させ,
その結果をアニメーションさせています.
移流方程式は微分方程式なのでそこで微分を使っています.
コンピュータは数学的に厳密に微分を計算できません.
極限を取れないからです.
そうかといって全く何もできないわけじゃなくて.
実は微分の定義にしたがって微分を近似して計算しています.
微分を近似すると実はただの引き算 (と割り算) になります.
さっきのページの上の方,
式 (2)-(4) を見てください.
こっちにも同じ式を書いておきますね.
\begin{align}
\frac{\partial u}{\partial x}
\approx
\frac{u(x+\Delta x)-u(x)}{\Delta x}.
\end{align}
\begin{align}
\frac{u_i^{n+1}-u_i^n}{\Delta t} + c \frac{u_i^n – u_{i-1}^n}{\Delta x}
=
0.
\end{align}
\begin{align}
u_i^{n+1}
=
u_i^n – c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n). \label{math_refuge00001}
\end{align}
見ての通り加減乗除の四則演算しかしていません.
特に上のページの (4) にあたる式 \eqref{math_refuge00001} に注目してください.
実はこの式にしたがって計算した結果をアニメーションでお見せしたのがさっきの動画です.
細かいことは追々やっていくとして何をしているか大雑把にいうと,
ベクトルなり数列なりの計算です.
また式 \eqref{math_refuge00001} を見てください.
右上にある添字 $n$ に注目すれば左辺は $n+1$ で右辺は $n$ です.
これは数列の漸化式とも思えますね?
数列も高校でやる内容です.
あなたが中学生ならちんぷんかんぷんでしょう.
でも実際に高校でやることがモノづくりなり何なりで
役に立っていることが感じてもらえればとりあえず OK です.
ちなみに移流方程式は上でもちょっと出てきた流体力学で出てくる微分方程式です.
もっと言うなら文字式の計算がゴリゴリ出てきていることも注意した方がいいでしょうか?
文字式できないとこの計算全く追えないですからね.
細かいことを言い出すときりがないので,
いったん今回はこのくらいにしておきましょう.
まとめると今回は次のポイントを掴んでもらえば十分です.
今回くらいのレベルでもベクトルで言うと
100 次元ベクトルとかそういうレベルの計算が必要です.
手計算ではやってられないのでプログラムを書いて計算して,
さらにプログラムを書いて図にしたりアニメーションにしたりしています.
このプログラミングについてはどこまで深掘りするかは未定です.
要望が多いなら別にきちんとやろうかとも思っています.
今回お見せしたのは変数が 2 つある偏微分方程式でちょっと難しいです.
変数が 2 つあるとややこしいので次回以降,
しばらく変数が 1 つしかない常微分方程式というのを見ていく予定です.
最後にアンケートをお願いしています.
改善に繋げていきたいのでぜひコメントをお願いします.
では次回をお楽しみに!