カテゴリー: プログラミング

  • プログラマと数学の接点, 入口: Python で数学を 番外

    参考にするためまずはメモ.
    どういうコンテンツ作っていこう?
    とりあえずは Sympy 頑張ろう.

  • 今後の中高数学勉強の指針/中高数学駆け込み寺 第 11 回

    はじめに

    こちらに 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 入門といった記事はあります.
    ただ他のサイトのコンテンツの方が網羅的で取り組みやすいのが正直なところです.
    まだそこまできちんと整備していないので.

    少し前まではドットインストールをお勧めしていたのですが,
    全部勉強するのに有料会員になる必要が出てきてしまいました.
    こちらも少しずつ情報を整備していく予定です.

    息抜きの数学コンテンツ

    漫画や小説からノンフィクション,
    数学者のエッセイまでいろいろ入れてあります.
    数学的に異常な難易度を誇る本も入れてあるので注意してください.
    実際に私が息抜きに気軽に適当に眺めている本で,
    基本的にどれもお勧めできるクオリティです.

    • [@KenjiFukaya3] 『数学者の視点』
      • Amazon へのリンク: http://tinyurl.com/zkrytkg
      • コメント: http://wp.me/p4PcgX-9V に書評を書いた.
        スーパー面白いのでとにかく買って読もう.
        学部 1 年のときにこの本を読んで数学への憧れを深めた.
        私の幾何への憧れもこの深谷賢治先生への憧れから来ている.
    • [@SuugakunoTanoshimi1, SuugakunoTanoshimi2, SuugakunoTanoshimi3] 『数学まなびはじめ』
    • [@KoujiShiga1] 『無限からの光芒―ポーランド学派の数学者たち』
      • Amazon へのリンク: http://tinyurl.com/zc4avpv
      • コメント: 文句なしで面白い.
        個人的に院で数学に行ったことに関して深い影響を与えている.
        関数解析が好きなら確実にはまる.
        作用素環の重鎮, 竹崎正道先生も絶賛しているレベルなので, とりあえず買って読んでおくべき.
    • [@NobukoUemura1] 『たかが数学, されど数学』
      • コメント: 山形大数学科の数学エッセイコンテストで入賞していた作品.
        理学部や数学科の HP 改訂でどこにあるのかわからなくなってしまった.
        数学エッセイの過去ページも見当たらない.
        人類の損失レベルの素敵な文章なのでどうにかしてほしい.
        そのうち山形大学の数学科に問い合わせたいと思っている.
    • [@KunihikoKodaira1] 『新・数学の学び方』
      • Amazon へのリンク: http://tinyurl.com/poysx3d
      • コメント: 学部 1 年のとき, もとの『数学の学び方』を読んだ.
        「こんな風にやるのか」と思って, 素直に実践していった.
        今思うと全然書かれている風にできていなかったが,
        それでも小平スタイルの勉強法に学部 1 年で触れられたのはよかったと思っている.
        願わくば中高生のときに知りたかった.
        そしてそんな気持ちがあったからこそ受験関係のプロジェクトをはじめた.
    • [@YasutakaIhara1] 『志学数学 -研究の諸段階 発表の工夫』
      • Amazon へのリンク: http://tinyurl.com/juubntl
      • コメント: 数学会の会誌で河東先生が「とりあえず買って読むべき」と書いたほどいい本.
        数学者の卵への思いやり溢れる穏やかな筆致で話が進んでいく.
        一般の人が「数学者はこんなことを考えながらこんなことをしているのか」という感じで読んでいっても十二分に楽しめるだろう.
    • [@Suurikagaku1] 『数学の道しるべ―研究者の道とは何か』
      • Amazon へのリンク: http://tinyurl.com/hyzotoz
      • 系統としては 『数学まなびはじめ』 [@SuugakunoTanoshimi1, @SuugakunoTanoshimi2, @SuugakunoTanoshimi3] と同じ.
        個人的には『数学まなびはじめ』の方が好きだが,
        こちらも十分楽しい.
        数学者が何を考えてどんなことをしているか知りたいならぜひ読んでみてほしい.
    • [@Suurikagaku2] 『物理の道しるべ―研究者の道とは何か』
    • [@HiroshiYuuki1] 『数学ガール』
      • Amazon へのリンク: http://tinyurl.com/glco42e
      • コメント: 中高生でも読みやすいと評判.
        ファンも多いし, 外国語への翻訳すらある.
        これを彼氏からプレゼントされて,
        算数・数学の勉強をはじめたという成人女性もいたりする.
        シリーズなので, 1 冊読んでみて合うようならいろいろ眺めるといいだろう.
        私はといえば, 「こんな学生生活なかったな?」という感じで,
        小説部分だけ眺めて数学部分をほとんど読み飛ばしていて,
        数学的内容はほとんど全く頭に入っていない.
    • [@MiyukiKawamura1] 『多面体の折紙 正多面体・準正多面体およびその双対』
    • [@OgataOzawa1] 『東大数理ビデオアーカイブス』
      • 東大ビデオアーカイブスへのリンク: http://tinyurl.com/j353cpx
      • コメント: 2015 時点で京大 RIMS 教授の小澤先生のアロハ装の真相について,
        当人の発言動画:3:40くらいからご自身で理由を語っている.
        東大数理ビデオアーカイブス http://www.ms.u-tokyo.ac.jp/video/ 自体面白い講義が揃っているので,
        興味がある向きはぜひいろいろ見てほしい.
    • [@KunioYasue2] 『量子の道草―方程式のある風景』
      • Amazon へのリンク: http://tinyurl.com/zq9t9zm
      • コメント: 方程式を絵画のように楽しんでみようという企画.
        目のつけどころが面白いし, 私もちょっとやってみたい.
        著者の妙に自慢気な文章スタイルが鼻につく人もいるだろう.
        私もたまにいらっとする.
    • [@HiroshiEzawa1] 『だれが原子をみたか』
      • Amazon へのリンク: http://tinyurl.com/zxe6p6u
      • コメント: 江沢節とでも言うべき異常なくらい力強い文章で原子を見ることを通じた物理での世界との向き合い方を考える本.
        楽しい.
    • [@YokoOgawa1] 『博士の愛した数式』
      • Amazon へのリンク: http://tinyurl.com/jrw84kl
      • 数学会の会誌の書評で「自分がやっている数学という営みが本当に素晴らしいものなのだと改めて気付かせてくれた」とあったくらい,
        数学者の心をも掴むよい本.
        超お勧め.
    • [@WilliamDumham1] 『微積分名作ギャラリ――ニュートンからルベーグまで』
      • Amazon へのリンク: http://tinyurl.com/hwcyuof
      • コメント: 中身はゴリゴリの数学ではあるが,
        有名だがあまり詳しくは扱わない色々な例・反例を紹介していたりするので,
        副読本として非常に勉強になる.
        Weierstrass の連続だがいたるところ微分できない関数の構成と
        証明が書かれていたりする.
        反例の構成や証明は大事なので,
        そういうところでも役に立つ楽しい本.
    • [@Hardy-Littlewood-Polya1] 『不等式』
      • Amazon へのリンク: http://tinyurl.com/huku9l4
      • コメント: 数学関係者向け.
        色々な不等式を扱った本.
        離散的な和から始まり, 積分不等式なども議論する.
        最近の論文にも引用されることがある程有名な本で,
        持っているだけでも幸せな気分になれる.
        必ずしも読み込む必要はないが,
        ちょっと変わった観点の数学入門という感じで数学を楽しみたい人にはいいだろう.
    • [@AiglerZiegler1] 『天書の証明』
    • [@YoshifumiTsuchimoto1] 『xのx乗の話』
      • Amazon へのリンク: http://tinyurl.com/hqf6q45
      • コメント: あっさりしているので何となく何が問題で
        どんな世界が展開されていくかを知るのには都合がよく,
        そういう読み方をするなら非常に面白い.
        もちろん, かっちり読むのにはつらいがそういう本ではない.
    • [@TatsuoKimura1] 『佐藤幹夫の数学』
      • Amazon へのリンク: http://tinyurl.com/hylx72j
      • コメント: 佐藤幹夫に関する色々な文章が載っている.
        単純な数学の話ばかりではなく,
        指導教官と (元) 学生とのほのぼのとした数学的な対話みたいな文章もあって楽しい.
    • [@KawazoeAi1] 『白と黒のとびら オートマトンと形式言語をめぐる冒険』
      • Amazon へのリンク: http://tinyurl.com/gsgk4tz
      • コメント: 小説だが文句なく面白い.
        オートマトンを迷路形式で解説している.
        まさにとりあえず買って読めというレベル.
    • [@ShoshichiKobayashi1] 『顔をなくした数学者–数学つれづれ』
      • Amazon へのリンク: http://tinyurl.com/z3m5dzk
      • コメント: 小林昭七先生の遺構となったエッセイ.
        完全版を読んでみたかった.
        これまた面白い [@SuugakunoTanoshimi1, SuugakunoTanoshimi2] (『数学まなびはじめ』) の昭七先生の記事の感想を書いたら弟の小林久志先生からご連絡を頂いて驚いたことがある.
        昭七先生が亡くなったのに合わせて昭七先生が執筆された文献を整理していて,
        検索したら『数学まなびはじめ』の私の昭七先生の感想記事を見つけたとかいう経緯だ.
        すぐに日本には戻れないので取り急ぎ記事のコピーを読みたいというのでお送りしたところ,
        とても喜んで頂けたのでほっこりした.
        何でも情報は出しておくものだと思った一件.
    • [@NaoyukiUchimura1] 『古都がはぐくむ現代数学: 京大数理解析研につどう人びと』
      • Amazon へのリンク: http://tinyurl.com/hh4c7hp
      • コメント: RIMS こと数理解析研究所に関する本.
        数学者の息吹を感じる.
        とにかく面白い.
    • [@KokuritsuTennmondai1] 『理科年表』
      • Amazon へのリンク: http://tinyurl.com/hahdqmq
      • コメント: 毎年更新するのもめんどいのでここでは平成 25 年度版を紹介しているが,
        新たに買うなら最新版を買おう.
        私は何となく持っていると格好よさそうという理由だけで高校 3 年のときに購入した記憶がある.
        パラパラと眺めているだけでも楽しい本だ.
        具体的なデータに親しみ,
        数値への感覚を持っておくのはとても大切なので,
        暇なときに何となく眺めてみる癖をつけてもいいくらい.
        楽しくデータを見て想像を膨らませてほしい.
    • [@HalTasaki6] 『数学:物理を学び楽しむために』
      • 文献へのリンク: http://tinyurl.com/ybwrkgw
      • コメント: 随時追記・改訂されている数学の教科書.
        例えば物理数学のシリーズも多くあるが,
        1 人の著者が責任を持ってきちんと統一したスタイルで書くべきだという信念に沿って執筆されている.
    • [@KentaroSato1] 『炭素文明論 「元素の王者」が歴史を動かす』
      • Amazon へのリンク: http://tinyurl.com/hnze4dc
      • コメント: タイトル通り有機化学に関する本だが抜群に面白い.
        こんなタイプの科学史の本が増えると,
        理工系の生徒・学生がもっと楽しく世界史や地理を勉強しやすくなる.
        私もこういうのを書いてみたい.
        死ぬ程時間かかるが.
        とりあえずは既にあるよい本はどんどん紹介したい.
    • [@TodaiKeisuu1] 『数理工学への誘い』
      • Amazon へのリンク: http://tinyurl.com/lohj6nw
      • コメント: 数学を工学的にどう使っていくかを説明している.
        いわゆるバリバリの物理系の話ではなく,
        「こんなところにも使われているのか」という系統の話題が載っている.
        「数学なんてどこで使うんだ」という人は読んでみるといいだろう.
    • [@KazuoHaga1] 『オリガミクス 幾何図形折り紙』『オリガミクス 紙を折ったら、数学が見えた』
      • Amazon へのリンク 1: http://tinyurl.com/hxx4r8b, Amazon へのリンク 2: http://tinyurl.com/zgy6t3q
      • 折り紙で数学するという不思議な毛色の本.
        この本ではないが, 折り紙と作図問題という研究もあり,
        動画などで紹介していきたいと思っている.
    • [@SimonSin1, SimonSin2] 『暗号解読 上下』
      • Amazon へのリンク 上: http://tinyurl.com/hgld6h9, Amazon へのリンク 下: http://tinyurl.com/gph68bg
      • コメント: 文句なしに面白い.
        ノンフィクションで暗号の歴史を追いかけながら,
        その背後にある数学にも迫っていく.
        これに限らずサイモン=シンの本はどれも面白い.
    • [@ScottPakin1] 『The Comprehensive \LaTeX Symbol List』
    • [@AitoAoyagi1] 『浜村渚の計算ノート』
      • Amazon へのリンク: http://tinyurl.com/j657jfg
      • コメント: 数学をネタにした推理小説.
        数学の扱いが悪いから数学者が反乱を起こした,
        というような話で数学サイドが悪者扱いになっているとも言える.
        いろいろな意味でときどき悲しくなるが,
        そうは思わず楽しんでいる人も多いようだ.
        数学をネタにした小説だし推薦しておく.
    • [@KeithDevlinGaryLorden1] 『数学で犯罪を解決する』
      • Amazon へのリンク: http://tinyurl.com/huummak
      • コメント: アメリカの刑事ドラマ Numbers で出てきた犯罪捜査に使われた数学を紹介している.
        単純な数学ばかりではなくセキュリティや暗号などの諸科学と関係の深い数学も紹介している.
        数学について深い解説があるというより数学と関わる広い世界を紹介しているような本.
        役に立つ数学に興味があるならお勧め.
    • [@YuukiOujou1] 『青の数学』
      • Amazon へのリンク: http://tinyurl.com/hagcm64
      • コメント: 数学が好きな高校生を集めて数学の決闘をしたりするという,
        ちょっと不思議なコミュニティを作っている天才数学者がいて,
        そこに参加している高校生達が描くストーリーを語る小説.
        最初はどんなもんかと思っていたらどんどん引き込まれた.
        数学それ自体に対する記述はほとんどない.
        むしろ数学に挑む人々の心のありようを描いた作品.
        決闘というと物騒だが, 相手を負かすことに力を注ぐよくあるタイプの人間もいれば,
        相手は関係なくただ数学と向き合う人間もいれば,
        「数学をこんなにも愛している人がたくさんいる」と思いながら決闘に挑む者もいる.
        数学の細かい話が出てくるわけではないので,
        数学が好きな人達の交流みたいなところに興味がある人には特にお勧め.
    • [@MasaeYasuda1] 『数学女子』
      • Amazon へのリンク: http://tinyurl.com/zyzpedu
      • コメント: 数学者が見ても楽しい漫画.
        作者は実際に数学科の学生だった人で,
        鹿児島大とその教官陣が本当にモデルになっているとのこと.
        飯高先生が愛してやまなすぎで雑誌「数学セミナー」で著者と対談したくらい.
        確かに「数学科あるある」濃度は極めて高く面白い.

    アンケートをお願いします

    今回もアンケートがあります.
    改善につなげるためぜひ回答をお願いします.

    この講座はいったんここで終了です.
    ここまでお疲れ様でした.
    細かい部分をあまり説明していませんし,
    けっこう大変な内容だったと思います.

    今後も数学や物理の学習に関する情報は配信していきます.
    ぜひお付き合いください.

    \clearpage

  • 漸化式から微分へ: 微分方程式のシミュレーションの観点から/中高数学駆け込み寺 第 8 回

    こちらに PDF があります.
    サイトでは見づらい方, コンテンツを手元に置いておきたい方は
    ダウンロードしてご覧ください.

    漸化式はルールの 1 つ

    前回数列の説明をしました.
    関数は数と数の対応ルールのことで,
    数列は関数の特別な形,
    つまり自然数と実数の対応ルールのことでした.

    で, 漸化式.
    漸化式は数列のルールの作り方の 1 つなのでした.
    ここから微分方程式から見た微分の発想につながっていきます.

    数列 $(a_n)$ の添字 $n$ の意味

    この添字, 普段は数列の $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}$ が曲者です.
    これを何とかしたい.
    もちろん一般的にはどうにもなりません.
    でも特殊な状況なら何とかなるかもしれません.
    数列と漸化式を基本に考えていても苦しいので,
    いったんそこから抜けることにしましょう.

    最初に漸化式は微分方程式を近似して出てくることを見ましたね?
    このミニ講座では微分方程式が本体です.
    何の指針もなく五里霧中を彷徨うよりも,
    大事な微分方程式と関係した漸化式を考える方が手が打ちやすいです.
    というわけで数列のレベルでぐちゃぐちゃ考えるのはこの辺にして,
    もう微分に行ってしまいましょう.

    今回はこの辺で終わりにします.
    お疲れ様でした.

    アンケートをお願いします

    今回もアンケートがあります.
    改善に役立てたいのでぜひご協力をお願いします.

    • https://goo.gl/forms/D0ZGnM7RADXoBJYH3

    ではまた次回をお楽しみに!

  • 素数夜曲とプログラミング B.2-B.4 Scheme, Python, Haskell, JavaScript 第 2 回

    素数夜曲とこれまでの内容

    素数夜曲

    これまでの内容

    B.2 名前と手続き

    リストの一番左は特別: そこには手続きを置く.
    命名と抽象化について割といいことが書いてある気がする.

    太字で強調してあるところがあって,
    何となく大事そうだからこちらにも書かせてもらおう.

    プログラミングは命令的知識を扱い, 数学は宣言的知識を扱う.

    命令的というのは「如何にして為すか」という意味で,
    宣言的というのは「何であるか」という意味.

    B.2.1 関数の定義

    素数夜曲では「函数」と書いているがめんどいので関数と書くことにする.
    いわゆるラムダの話から入る.
    次の関数は $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!"
    

    Python の lambda

    $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

    Haskell での lambda は次の通り.

    
    mysquareInt = \x -> x ^ 2
    mysquareDouble = \x -> x ** 2
    
    main = do
      print $ mysquareInt 2
      print $ mysquareDouble 2.0
    

    RESULT

    4
    4.0
    

    JavaScript の lambda

    一般的な 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
    

    B.2.2 アルファ変換

    いわゆる束縛変数の文字を変えるというやつ.

    B.3 ラムダ算法

    Lambda calculus だし,
    たぶんいわゆるラムダ計算なのだろう.
    プログラムじたいはないのでばっさり省略.

    B.3.3 イータ変換

    省略.

    B.3.4 ラムダ項の定義

    略.

    B.3.5 コンビネータ

    略.

    B.3.6 簡約の戦略

    略.

    B.4 特殊形式 (special form)

    対象に名前をつけてプログラム上のどこからでも使えるようにすることを,
    トップレベル定義 (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
    

    アルファ変換は局所変数を別の文字 (列) に置き換えてもいいとか,
    とりあえずそういう雑な理解をしている.
    計算機科学的にもっと面倒なところまでカバーしている可能性があり,
    そこまで調べていないからいまの理解の雑さは割と真剣に警戒している.

    B.4.1 define

    用語の定義があるから適当に書いておこう.
    対象を挿入する場所 (スロット) は本文と同じくを角括弧 (<>) で書きたいが,
    いろいろな事情から隅括弧 ([]) で書くことにする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 での関数定義を書いておこう.

    Python での変数・関数定義

    まずは変数定義.

    
    a = 1
    print(a)
    

    RESULT

    1
    

    数を 2 乗する関数を定義する.

    
    def f(x):
        return x ** 2
    
    print(f(2))
    print(f(3))
    

    RESULT

    4
    9
    

    Haskell での変数・関数定義

    私が理解している限り, いわゆる変数はない.
    定数関数があってそれが定数の役割をしてくれるという理解.

    
    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
    

    JavaScript の変数・関数定義

    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: procedure
    • exp: expression
    • num: number
    • predi: predicate
    • func: function
    • var: variable
    • int: integer
    • consq: consequent
    • args: argument
    • val: value
    • str: string
    • altna: alternative
    • char: character
    • obj: object
    • lst: list
    • stm: 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) と呼ぶ.

    B.4.2 lambda

    
    定義方法: (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 strictlet, const を使うのがいいのだろうか.
    識者のご意見求む.

    begin

    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 で似たようなのがあるだろうか?

    今回はこの辺で終わろう.


    1. 実体参照にすればいいのだが,
      今の私の腕ではスロットだけを自動で変換するスクリプトを作れないので諦めた.
      Markdown の引用の「>」や TeX 中の不等式を無視して,
      スロットだけを狙い打って実体参照に書き換えるコードを書く腕がない. 
  • 素数夜曲とプログラミング B-B.1.2 Scheme, Python, Haskell, JavaScript 第 1 回

    はじめに: なぜ素数夜曲か

    数学・物理学習とプログラミングと絡められないかと試行錯誤している.

    数学からの取っつきやすさ,
    物理からの取っつきやすさ,
    プログラミングからの取っつきやすさをそれぞれ考えないといけない.
    ある程度の体系性もほしい.

    数学や物理としても面白い内容にしたいし,
    プログラミングから見ても意味がある内容にしたい.

    どうしたものかとずっと思っていたのだが,
    素数という数学のキラーコンテンツとプログラミングを絡めた素数夜曲があった.

    微妙なところも多いとは思いつつ,
    メインエディタとして 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
    

    B.1 評価値を得る

    とりあえず Hello, World! で.

    
    (write "Hello, Scheme!")
    

    RESULT

    "Hello, Scheme!"
    

    B.1.2 四則計算

    基本的な書き方

    加減乗除は次の通り.
    前置記法なのがポイント.
    割り算で分数にしてくれるのは割とポイント高い気がする.

    
    (write (+ 5 3))
    (newline)
    (write (- 5 3))
    (newline)
    (write (* 5 3))
    (newline)
    (write (/ 5 3))
    

    RESULT

    8
    2
    15
    5/3
    

    Python

    
    print(5 + 3)
    print(-5 + 3)
    print(5 * 3)
    print(5 / 3)
    

    RESULT

    8
    -2
    15
    1.6666666666666667
    

    Haskell

    まだ (?) 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
    

    JavaScript

    
    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

    Python だと普通に書く.

    
    a = 2 / 3 + 5 * 7 - 11
    print(a)
    

    RESULT

    24.666666666666664
    

    Haskell

    Python っぽくも Scheme っぽくも書ける.
    演算子を括弧でくくると Scheme のように演算子を前に置けるから.

    
    main = do
      print $ 2 / 3 + 5 * 7 - 11
      print $ (-) ((+) ((/) 2 3) ((*) 5 7)) 11
    

    RESULT

    24.666666666666664
    24.666666666666664
    

    JavaScript

    
    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 版は次の通り.
    Python は math をインポートしないといけないのがめんどい.
    Python 3.5.1 を使っているが関数合成はないのだろうか.
    ちょっと調べたところでは標準ではなさそうだった.

    
    import math
    print(math.sqrt(2))
    print(math.sqrt(math.sqrt(4)))
    

    RESULT

    1.4142135623730951
    1.4142135623730951
    

    Haskell

    Haskell は関数合成 (.) がある.

    
    main = do
      print $ sqrt 2
      print $ sqrt $ sqrt 4
      print $ (sqrt . sqrt) 4
    

    RESULT

    1.4142135623730951
    1.4142135623730951
    1.4142135623730951
    

    JavaScript

    
    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

    Python の場合はやはり math モジュールにある.

    
    import math
    print(math.e)
    

    RESULT

    2.718281828459045
    

    NumPy にもある.

    
    import numpy as np
    print(np.e)
    

    RESULT

    2.718281828459045
    

    Haskell 版

    $e$ 自体は定義されていないので Scheme と同じく指数関数から近似値を出す.

    
    main = do
      print $ exp 1
    

    RESULT

    2.718281828459045
    

    JavaScript 版

    Math.exp() が自然対数の底による指数関数で,
    Math.pow() は一般の底と指数の関数.

    
    console.log(Math.exp(1));
    console.log(Math.pow(2, 3));
    

    RESULT

    2.718281828459045
    8
    

    2 項演算子をいくつか

    まずは等号 = を.

    
    (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
    

    Python

    等しくないことの判定は != だ.

    
    print(1 == 2)
    print(1 == 1)
    print(1 != 1)
    print(1 < 2)
    

    RESULT

    False
    True
    False
    True
    

    こちらは = だと代入になってしまうので == と重ねる.

    Haskell

    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
    

    Haskell 追記

    またコメントを頂いてしまった.
    これ.

    
    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 もきちんと全部読んだわけではないし,
    この雑な理解でいいとも思わないがとりあえずこれで掴まえられるところもある.
    「適当にしか理解していない (そして半端な理解だと割とどうでもいいところで無駄にはまる)」ことさえ理解しておけば,
    あとでいくらでも調節は効くと数学でいつも経験することをいつも通り信じて追記を終える.

    JavaScript

    比較には ===== がある.
    基本的に === でやるのが安全.

    
    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 の仕様は改定を重ねられているので、どの版であるかによって挙動が異なる部分があります。 Scheme の仕様はその正式名称を略して RnRS と呼ばれていて、 n の部分に何回目の改定であるかの番号が入ります。 最新は R7RS です。
    「大文字小文字を区別しない」というのは R5RS までの仕様です。 (ただ、 R5RS 処理系を名乗る処理系でも区別する処理系や、切り替えられる処理系はあります。) R6RS 以降では区別するのがデフォルトです。
    有理数や無限長整数などは R6RS 以外ではオプショナルな仕様です。 (R6RS で必須になった後、 R7RS では再びオプショナルになりました。) 処理系によっては桁あふれに配慮が必要な場合もあります。 有理数や無限長整数は Guile や Gauche では問題なく使えますので気にする必要はありませんが、その他に色々なところで仕様では未規定の部分があるので処理系を乗り換える機会にはつまらないことで躓くかもしれません。

    思っていたよりも面倒そうな話だった.
    ただ「ハマる可能性があるかもしれない」
    「ハマるところはこの辺」というのが頭の片隅にあるかないかだけでも
    全く違うことはこれまでの経験でよくわかっている.
    こんなありがたいことはない.
    こういうコメントを頂けるのが記事にしておく醍醐味だ.

    地道に続けよう.

  • 数列とは何か? 微分方程式のシミュレーションの観点から/中高数学駆け込み寺 第 7 回

    こちらに 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 つ前と関係あるようにしてもいいです.
    この一般的なやつがいわゆる漸化式です.

    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 回.

    一番単純でしかも実際に使われる微分方程式としてまずは 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}

    これに沿って計算したのがいわゆるオイラー法.
    次のセルではこれをコードに落としている.

    1 階の方程式のプログラム

    
    %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 階の常微分方程式: 単振動

    こちらは単振動.

    次に 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

    ルンゲ-クッタ短評

    今度の一致具合はなかなかよさそう.
    あくまで見た目の感じではあるけれども.

    良くなったのは一番下の図, 相空間軌道だ: ちゃんと閉じてくれた.

  • ベクトルとは何か? 微分方程式のシミュレーションの観点から/中高数学駆け込み寺 第 5 回

    こちらに 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')
    
  • 経済や生物で使う微分方程式/中高数学駆け込み寺 第 3 回

    こちらに 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


    1. 本当は $x(t)$ は整数なんですがここではいったん実数だと思っておきます.
      この辺の処理は物理や化学でもよく出てきます.
      実はこの近似というか合理化もきちんとやっておかないといけません.
      ここではとてもやりきれませんけど.
      本当は時間についてもいろいろ議論があります. 
    2. もっと広く力学系と呼ばれる分野にまで広がります.
      力学系まで行くととんでもなく難しくなるのでとてもここで紹介はできません.
      でも幾何学とも関係してきたり整数論とも関係してきたり,
      射程範囲は広いです. 
  • 厳密解との比較: 放射性物質の崩壊/中高数学駆け込み寺 第 2 回

    こちらに 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 等分なのだ.

    ここで刻みが荒いといい近似にならないことを大雑把に調べておこう.

    2 等分

    
    
    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

    4 等分

    何となくそれっぽい形にはなっている.
    もちろん精度は大したことない.

    
    
    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

    10 等分

    あまり適当なことを言うのも良くないが, けっこう滑らかに見える.
    精度はまだまだ.

    
    
    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

    50 等分

    まだ差が目で見える.

    
    
    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

    1000 等分

    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 つしかない常微分方程式というのを見ていく予定です.

    最後にアンケートをお願いしています.
    改善に繋げていきたいのでぜひコメントをお願いします.

    では次回をお楽しみに!