計算数理A演習第1・2回

第1・2回(第1週)


はじめに

計算数理A演習では、プログラミング言語としてPythonを用います。 計算数学演習ではコンパイルを必要とする言語であるC言語を用いていましたが、 本講義ではインタプリタ言語と呼ばれるPythonを用います。

Pythonは近年もてはやされているディープラーニングによってかなり有名になりました。 その理由の一つに、言語自体の設計の新しさと手軽さやライブラリの豊富さがあります。

プログラミング言語は慣れてしまえば新しい言語を学ぶことは難しくないですが、初級者である皆さんに二つの言語を学んでもらうことはリスクでもあると思います。

しかし、Pythonの手軽さは新しい言語を学ぶコストよりも手軽に出来ると思いますし、C言語の学習で学んだfor文・if文などの考え方は全く変わりません。 標準語が関西弁に変わったような感覚で書き換えていける程度の差しかありません。 ですのでアルゴリズムについてわからないところがあれば、 計算数学演習のHPはこちらにありますので、参考にしてください。


今週の演習

今日は以下の内容で演習を行います。

  • 本ターム内に、この演習で行うこと
  • 演習のホームページの紹介と見方
  • Pythonの導入

本ターム内に、この演習で行うこと

藤本先生が担当される計算数理Aの講義には毎週必ず出席してください。 この演習の内容は計算数理Aの講義内容に強く依存し、講義内容をよく理解している事を前提条件としています。 演習のみ登録されている人も、できれば講義にも出てください。

演習の当面の目標は、コンピュータを使ってある種の常微分方程式・偏微分方程式を解くことです。 ただし、その為には

  • プログラミング言語の習得
  • 計算結果の可視化技術の習得
  • 数値計算法の習得

といったことが必要となります。

本演習では、Pythonを用いますが、基本的な数値計算の方法は計算数学演習で学んだことがベースになります。

演習序盤(次回から2~3回程度)はPythonの使い方と、グラフィックスライブラリの使い方の演習を行います(予定)。

次に授業で紹介される常微分方程式を数値計算を用いて解いてもらいます。 最終的には、拡散方程式(熱方程式とも呼ばれる)や波動方程式といった偏微分方程式を数値計算で解いてもらい、結果を可視化してもらいます。

演習では、数値計算法の習得も目標の一つではありますが、数値計算の必要性や、それを道具として使うという姿勢について、ある程度理解してもらえればと思います。


演習のホームページの紹介と見方

演習のホームページのトップに「お知らせ」という項目があります。 そこには、レポートの出題状況やその他大事なお知らせを載せますので、毎週必ず見てください。



Python入門

Hello, World!!

C言語でも初めに行うHello, World!!の書き出し方法です。 エディタに次を書き込んで、main.pyとして保存してください。

print("Hello, World!!")

そして、

python3 main.py

とターミナルで実行してください。

※ Pythonには現在2系と3系と呼ばれる異なるバージョンが広まっています。最近では標準で最新の3系がインストールされますが、2系のバージョンを使用している場合プログラムが回らないこともあります。もしも上記方法でプログラムが回らない場合には、ターミナル上で
python -V
と入力してバージョン確認を行い、3系のバージョンへ更新してください。

比較のため、C言語で同じものを書くために必要なプログラムを載せます。

#include <stdio.h>

int main() {
    printf("Hello, World!!\n");
    return 0;
}

このプログラムに加えて、コンパイルと実行の二つの処理がさらに必要になります。

このようにコードとして書く量を大幅に減らすことができ、結果的にミスを少なく出来るのがこの言語の大きな特徴になります。

さて、この簡単なプログラムからPythonとC言語の二つの言語に見られる違いを述べておきたいと思います。

    • Pythonでは文末に ;(セミコロン)を書いてはいけない。
    • main関数が存在しない。

Pythonにはmain関数がないのですが、一応それに対応する考えはありますがここでは触れません。


Spyderを使う

 さらに大きな武器を使ってプログラムの学習を進めていきたいと思います。 これまではエディタでプログラムを書き、ターミナルでコンパイルして実行してファイルにデータを保存して、gnuplotでグラフを書きました。  本演習では、これらのツールがひとまとめになっているソフト(統合開発環境)であるSpyderを使います。言語としてはPythonは言語の書き方としてC言語にかなり近いです。

※ 2019年度の演習ではJupyter Notebookを使いました。こちらはMathematicaと似て、ざっくり言うとエディタとターミナルが入り混じったような形式でした。SpyderとJupyter Notebookは、どちらも一長一短です。本演習では、基本的に受講者のプログラミング経験は計算数学演習でC言語を学習しただけ、として想定していますので、エディタが分かれていて感覚的にプログラミング感覚が近いSpyderを利用することにします。ただし、開発環境の違いは、本演習の本来の目的には影響しません。どちらの形式で進めてもらっても構いませんが、

レポートなどの形式は「〇〇.py」の形のPythonプログラム形式(〇〇.ipynbではない)とします。

それでは、Spyderを実行しましょう。

macOS、Linuxの場合

まずターミナルを開いて、

spyder &

と実行してください。

Windowsの場合

検索欄に”anaconda prompt”と入力し、検索結果からアナコンダプロンプトを開きます
(下図参照)。

このアナコンダプロンプトがターミナルの役目を果たします。

もしspyder &コマンドでエラーになる場合は” &”を取り除き、spyderとだけ記述し実行してください。

しばらくすると次の様な画面が現れるはずです。(自動更新に関するアラートはOKしておいてください)

左側がエディタ、右側上部がヘルプで、右側下部がiPythonコンソール(以降、コンソールと呼ぶ)です。コンソールは、とりあえずPython用のターミナルのようなものと覚えておいてください。このコンソール自体でPythonのコマンドを実行できますし、エディタで記述したプログラムファイルを読みこむこともできます。 コンソール下部のタブを選択することで、過去にコンソールに打ち込んだ記録(ログ)を参照することもできます。同様に、ヘルプ下部のタブを選択することで、定義されている変数一覧である変数エクスプローラーや、ファイルエクスプローラーを表示することができます。

Spyderのレイアウトは自由に変更可能です(デモ)。

プログラムの実行

それでは早速使ってみましょう。 デフォルトでは、現在いるディレクトリはホームディレクトリ(/fs/home/アカウント名)です。このままここで作業すると、ホームディレクトリ直下がごちゃごちゃすることになりますので、適当にディレクトリを作成(右クリック→新規→フォルダ)し、移動(ダブルクリック)してください。これでコンソール上でも移動したことになります。

試しに先ほど作成したmain.pyを実行してみましょう。 main.pyを作業ディレクトリに持ってきて、ファイルエクスプローラー上でダブルクリックしてください。すると、左側のエディタにmain.pyが表示されます。左上にある緑の三角(下図参照、バージョンによって違うかも)をクリックすると、main,pyが実行され、コンソール上に「Hello, World!!」が表示されます。

ファイルを作成するときは、左上の「ファイル」→「新規ファイル」で作成できます。この段階ではファイルは保存されていないので、「ファイル」→「保存」や「形式を指定して保存」などでファイル名をつけて保存してください。

コンソールの使い方

「Hello, World!!」と表示させたいだけなのであれば、いちいちファイルを使わなくても、interactiveにコンソール上で実行しても良いでしょう(直接コンソール上でprint("Hello, World!")と打ってreturn)。他にも、ちょっとした作業であれば、コンソールで事足ります。 ただし後々、行った作業をまとめて1つのファイルにする場合には、コンソール上で行なった作業も追加することを忘れないように。

コンソール上で複数行をまとめて実行したい場合は「ctrl + return」で改行できます。

コンソールとエディタの使い分けについて

Pythonのプログラムは、「コンソールでの直接入力」も、「エディタでの入力・実行」も同じ結果になります。 今後課題やレポートなどでプログラムを作成することが増えますが、どちらを使っても構いません。例えば、全てコンソールを使って作成しても良いですし、課題ごとに〇〇.pyファイルを作っても良いです。 ただしレポートでは、1個の完結した〇〇.pyファイルとして提出してください。したがって、全てコンソールで作成する場合は、必要な箇所を〇〇.pyファイルにコピーさせる必要があります。このとき必ず全ての変数が意図通りに定義できているか確認してください。提出されたファイルを実行したときに、変数の定義部分がなくて正しく実行できない場合、点数をつけることはできません。また、これ以降登場するnumpyやmatplotlibなどのモジュールの定義、importも忘れないようにしてください。

※ 以下では複数行に渡るコードも書かれていますが、コンソールでは基本的にブロック文以外では複数行で書くことができず、1行ごとの実行になります。


コメント文

Pythonでは、# 以降の部分はコメント文になります。
C言語における // と同じです。

print("これは表示される")
#print("これは表示されない")
print("これは表示されるけど、コメントが後ろにある") # ここはコメント

変数

データ型と初期化

C言語との大きな違いが変数にも表れます。

  • 変数の型の”宣言が”存在しない

C言語ではint型やfloat (double)型の変数の宣言について学びました。例えばC言語では

float x;

のように、変数を使う前に宣言をする必要がありました。

Pythonではこのような宣言は不要です。
とはいえ、C言語のときと同様に変数の型は存在し、データ型と呼ばれています。 また、宣言は必要なく初期化を以下のように行います。

x = 0.0

0.0と書けば小数点があるのでコンピュータが勝手にfloat型と認識してくれます。int型にしたければ、

x = 0

というように整数だけで初期化を行えば良いのです。

小ネタ:変数に代入されている値を確認

初期設定では右上にある区画。下のタブが”ヘルプ”になっていると思いますが、”変数エクスプローラー”のタブを選択すると、変数に代入されている値を見ることができます。非常に便利です。ちなみにここから変数を選択して、変数の値を変えることができます。

変数の出力

数値計算では、計算した変数の値を出力する必要があります。 C言語ではint型やfloat型に合わせて、%dや%lfなどと書く必要がありましたが、次のプログラムをCellに書き込んでRunを実行してください。

a = 1
x = 2.1
print(a, x)

変数がそのまま出力されたことがわかると思います。

小ネタ:print文で変数の前後の空白を消すには

例えば

a=1
print("a=", a)

とすると、

a= 1

となり、1の前に空白が入ります。場合によってはこのような空白が入って欲しくないこともあるでしょう。そのときには、

a=1
print("a=", a, sep="")

とすれば、空白がなくなり、

a=1

と表示されます。これは文字列と文字列の区切り(separater)を””で括られた文字列で置き換えることを意味しています。つまり、sep="\t"とすればタブで区切ることもできます。 今後、ケースバイケースで使い分けてみてください。

変数の入力

C言語ではscanf関数を用いましたが、以下のようにします。

s = input()

input関数は、どのようなものも入力できるため全ての場合に対応できる文字列として読み込まれます。 そのため、int型もしくはfloat型として読み込むには、以下のように型を指定して変換します。

s = input()
# int型
s = int(s)
# float型
s = float(s)
# 直接読み込む
s = int(input())
s = float(input())

Pythonでは初期化の際に型が決まるので、数値なのか文字列なのかを先に指定する必要はありませんが、何の指示もないとその後不便です。 scanfには文字が書けなかったので、printf関数が必要でしたが、Pythonでは以下のように出来ます。

s = input("sに数値を入力してください。s = ")

四則演算+α

Pythonでも四則演算はC言語と同様ですが、冪乗も演算子が用意されています。

演算 数学的表現 Python C言語
a + b a + b a + b
a - b a – b a – b
a \times b a * b a * b
a \div b a / b a / b
剰余 a \mod b a % b a % b
冪乗 a ^ b a**b または pow(a, b) pow(a, b) (要 数学関数)

a^bは冪乗演算子ではなく、bit演算子のXORに対応するので、冪乗と思って使うと異なる結果になるので注意!

また、次の代入演算子はPythonでも使うことができます。

代入演算子 Pythonでの意味
a += b a = a + b
a -= b a = a – b
a *= b a = a* b
a /= b a /= b
a %= b a = a % b
a **= b a = a**b

ところが、a++a--といった、C言語のfor文などでお馴染みだった表現は使えませんので注意です。

リスト

計算数学演習では、簡単にですが配列について学びました。 Pythonではリストと呼ばれ、以下のように定義します。 ここでも変数の宣言がないため、初期化を伴います。

X = [1, 3, 5, 7]

配列の要素へのアクセス、及び値の書き換えはC言語と同様に行えます。このとき、配列の番号は0から始まります。 3を出力したければ、

print(X[1])

5を書き換えて0にしたければ、

X[2] = 0

と出来ます。ただしこのやり方では100個の要素を持つようなリストを作ることは現実的ではありません。メインにはこの後に学ぶarray型というものを利用します。

小ネタ:履歴を確認

コンソールの下の方の”iPythonコンソール”と書かれているタブの隣に、”ヒストリ”のタブがあります。
その日やったことをまとめて、.pyファイルに保存するのに便利です。


制御文

インデント(字下げ)の注意点

C言語での制御文(if文、for文など)は、

for (i=0; i<10; i++) {
    ブロック文
}

のように { } で囲っていました。そのため、{ } 内のインデントが少々ぐちゃぐちゃでも、実行結果には影響を与えませんでした。 以下で学ぶPythonの制御文では、{ } の代わりにインデントで制御文で制御されるブロック文を判別します。そのため、インデントが正しく揃ってない場合、エラーや予期せぬ挙動をします。今一度インデントについて確認しましょう。 Pythonでのインデントは、スペース4個分が標準的です。Spyderでは、タブ(Tab)を押すとスペースが4個分の倍数の位置に止まりますので、タブを使うと良いでしょう。

ともかく、習うより慣れろ、ということで以下の例を見てみましょう。


if文

次にif文について学びます。条件式はifと : の間に挟まれ、制御文はif文のあとインデントされている行が全てif文で制御される範囲になります。

早速以下のサンプルコードを実行してみましょう。さらに、a=1に変えて実行し直してみてください。

a = 0
if a == 0:
    print(1, a)
    a = a + 2
    print(2, a)

a = a + 10
print(3, a)

if文の制御から外すには、if文の初めのところとインデントを揃える必要があります。

それでは、if文による分岐の大事な要素であるif-else文についてです。 先のサンプルコードに少し変更を加えるとif-else文が完成します。

先と同様に最初のa=0をa=1に変えてみてください。

a = 0
if a == 0:
    print(1, a)
    a = a + 2
    print(2, a)

else:
    a = -11
    print(3, a)

a = a + 10
print(4, a)

それでは、条件が複数ある場合についてです。 C言語ではif~else if~ else文でしたが、二つ目以降の条件がelifと省略されます。

a = 0
if a == 0:
    print(1, a)
    a = a + 2
    print(2, a)

elif a == 1:
    print(3, a)

else:
    a = -11
    print(4, a)

a = a + 10
print(5, a)

比較演算子と複数条件

以下の表はC言語と全く同じです。

演算子 説明 記述例 記述の意味
> 大きい if (a>10) aが10より大なら
>= 大きいか等しい if (a>=10) aが10以上なら
< 小さい if (a<10) aが10より小なら
<= 小さいか等しい if (a<=10) aが10以下なら
== 等しい if (a==10) aが10なら
!= 等しくない if (a!=10) aが10でないなら

しかし、複数の条件を条件式に入れる場合の論理演算子は以下のようになります。

演算子 対応するC言語演算子 記述例 記述の意味
and && if a==0 and b==1: aが0かつbが1
or || if a==0 or b==1: aが0またはbが1
not ! if not a==1: aが0ではない

以下の課題は計算数学演習でもやった課題。
アルゴリズムがわからなかったら過去のCで書いたプログラムを参照しましょう。

課題1

2つの実数 s, t を入力すると、その差の絶対値を表示するプログラムを作成せよ (if 文を使うこと)。

課題2

整数 m が入力されると、それが奇数か偶数かを判定するプログラムを作成せよ。

(注:剰余を求める演算子 % を利用. これは整数型変数m, n に対し、m%n は m を n で割ったときの余りを与える。この演算子は整数型変数にのみ利用可能。)

課題3

2つの実数 s, t を入力すると, s > t, s = t, s < t か判定して表示するプログラムを作成せよ。

課題4

3つの実数, a, b, c を入力すると, x に関する2次方程式, ax2 + bx + c = 0 が実数解を持つか判定し, 表示するプログラムを作成せよ。

課題5

入力された整数が、3 の倍数か 7 の倍数であれば、その数値を表示するプログラムを作成せよ。( or を使う。)

課題6

“3 の倍数か 7 の倍数”である三桁の自然数の入力を促し、入力された数値が 条件を満たしているかどうか画面に出力するプログラムを作成せよ。(if文は1回しか使えないものとする))


繰り返し文

for文

次に重要なfor文についてです。これはC言語とかなり異なります。 C言語のfor文では初期化式・条件式・反復式の三つがありましたが、Pythonのfor文は、

for 変数 in リストまたはrange型:
    文

のような構造になっています。例えば、0から9まで出力するには以下のような方法があります。

for i in [0,1,2,3,4,5,6,7,8,9]:
print(i)

for i in range(10): print(i)

2つ目の方法でrangeという関数を使っていますが、このrangeはいくつか使い方があります。

1. 繰り返し回数のみ指定する場合(i = 0 … N-1)

N = 100
for i in range(N):
    print(i)

このように1つだけ数字Nを入れると、0からN-1まで、繰り返していきます。

2. 開始と終了繰り返し回数のみ指定する場合

N = 100
for i in range(1, N, 3):
    print(i)

これを実行すると1, 4, 7, …, 97という数字が出力されるようになります。rangeの1つ目の引数が開始点、2つ目の引数が終了点+1、3つ目の引数が増加量を表します。C言語で考えると

int i;
int N = 100;
for (i=1; i<N; i+=3)
{
    printf("%d\n", i);
}

という風になります。もし増加量が1で良ければ最後に数値を入れる必要がなく、以下の二つのfor文は同じものになります。

N = 100
for i in range(3, N):
    print(i)

for i in range(3, N, 1):
    print(i)

while文

while文はインデントで制御文の範囲を指定する以外は、C言語に近いです。

例えば、

N = 100
i = 0
while i < 20:
    i = i + 1
    print(i)

小ネタ:print文の最後の改行を消すには

例えば

for i in range(5):
    print(i)

とすると、

0
1
2
3
4

となり、出力ごとに改行が入ります。

for i in range(5):
    print(i, end=", ")

とすれば、改行がなくなり、

0, 1, 2, 3, 4, 

と表示されます。同様に、

for i in range(5):
    print(i, end="")

とすれば、

01234

と連続して表示されます。

これは通常改行文字(”\n”)である文字列の終端(end)の処理を””で括られた文字列で置き換えることを意味しています。もちろん、上で紹介したsepと組み合わせて変えることもできます。

for i in range(5):
    print(i,"、の次の", sep="", end="")
print(i+1)

例題

C言語でも取り組んだ和と積を求めるプログラムは以下のようになります。

wa = 0
seki = 1
for i in range(1,5):
    wa = wa + i # wa += i でも可
    seki = seki * i # seki *= i でも可

print("和=", wa)
print("積=", seki)

課題7

和と積のプログラムを改良して, 整数 n を入力すると 1 ~ n までの和と積を表示するプログラムを作成せよ. (注意:あまり大きな n を入力すると、特に積の値がコンピューターで扱える範囲を超えてしまうので、正確な値がでてこなくなります. それほど大きくない n で試してください.)

課題8

実数 a, 自然数 n を入力すると an を計算し表示するプログラムを作成せよ. (注意:あまり大きな n を入力すると、特に積の値がコンピューターで扱える範囲を超えてしまうので、正確な値がでてこなくなります. それほど大きくない n で試してください.)

課題9

実数 a, 整数 n を入力すると a + a2 + a3 + … + an を計算し表示するプログラムを作成せよ.

課題10

自然数 n を入力すると n 以下の 3 の倍数と 7 の倍数を表示するプログラムを作成せよ.(for 文と if 文を使う. )

課題11

自然数 n を入力すると n が素数がどうか判定するプログラムを作成せよ.

以下、難しめの問題

課題12

三桁以下の(10進数の)自然数を入力したら、その数を二進数で表示するプログラムを作成せよ(for文を使って) 例: 三桁以下の自然数を入力してください:231 その数を二進数であらわすと11100111です。

課題13

7のn乗をAnとあらわすとする。下三桁が(943)となる最小のAnのnを求めよ。 そのようなnが存在しない場合は、そのように出力せよ。
ヒント:条件をみたすnが存在しない場合は、943を含まない集合で閉じている(循環している)ことになる。つまり、少なくともn=〇までに943がでなければ、条件をみたすnが存在しないと言って良い。

C言語では、intで扱える整数は、-2147483648から2147483647であったが、Python3では上限下限はなくなった。ただし「メモリが使える範囲で」という上限付きなので、プログラムが利用できるメモリの範囲を超えてしまうと、プログラムが異常終了するので注意。

課題14

整数Mのn乗をMnとあらわすとする。三桁以下の整数Mの入力を促し下二桁が(16)となる最小のMnのnを求めよ。そのようなnが存在しない場合は、そのように出力せよ。


関数

関数を利用者が自ら定義して利用することは、どのプログラミング言語でも必要な基礎的な技術です。早速Pythonにおける関数がどう定義されるかを見てみましょう。 まず、最も基礎的な入力された変数を書き出す関数です。

def myfunc(a):
    print(a)

returnがないのでC言語のvoid型の関数になります。しかしどこにもvoid型などと宣言するところがありません。 これは、変数の宣言がないのと同様に実際にreturnされる値が書き込まれるまでPythonでは関数が返す変数の型は定義されないので、変数と同様に必要がありません。 そして、:(コロン)で改行した後にはif文やfor文と同じように、インデントされています。これは理由も同じで、インデントされた部分だけが関数の中で処理される文と認識されるためです。 次に二つの和を返す関数を考えます。

def my_sum(a, b):
    c = a + b
    return c

a = 1
b = 2
ans1 = my_sum(a, b)
ans2 = my_sum(3, -5)
ans3 = my_sum(2.0, 4.1)
print(ans1, ans2, ans3)

このようにするとaとbの値が入力されると、その二つの変数の和をとって返してくれる訳ですが、引数にも型の定義がないため、実際に入力される変数の型をそのまま利用して和を取ることになります。そのため、ans1のような使い方に加えて、ans3のようにfloat型の値を直接入力することが出来ます。

注意点としては、C言語と違い先に関数を宣言しておくプロトタイプ宣言ができません。基本的に宣言という概念がないので当然ではあるのですが、常に利用するよりも前に関数を定義しておく必要があるので注意してください。


さて、ここまでに駆け足に基礎的な部分を学習してきました。

これだけではC言語でも同様にプログラミングができ、Pythonを学ぶありがたみを感じるには機能が少ないと感じるでしょう。 Pythonの強みは、その豊富なモジュールにあります (もちろん他にもありますが…)。 本講義で主に用いる、NumPyとMatplotlibというモジュールに使用方法を学び、実際に問題を解いてPythonに慣れていきたいと思います。

NumPy:数学関数、ベクトルなどの利用

C言語では様々な数学関数を用いるためにmath.hをインクルードしましたが、 それに対応するライブラリ(モジュール:後述)は、numpyとなり以下のように読み込みます。

import numpy as np

C言語で#includeに対応するものがimportになります。後ろにas npとありますが、これは、Numpyに含まれる関数を使う際に、「毎回numpy.〇〇と打つのが面倒なので、代わりにnp.〇〇を使う」という宣言です。numpyに関してはもはやセットでas npを使うというのが世界的にも慣例になっています。他のところから引っ張ってくる色々な関数や定数の集まりのことを、他の言語ではライブラリと呼ぶことがが多いですが、Pythonでは一般的にライブラリではなくモジュールと呼びます。特にNumPyはNumericalのNumから名前の由来のある数値計算用のモジュールで、これがあればかなりのことが出来ます。

数学関数・定数の利用

数学関数・定数などは以下のようにして使えます。

np.sin(x)
np.exp(-2*t)
np.pi

上記のxやtなどの変数は、1個のスカラーでも良いですし、リストでも構いません。 np.piはπを表します。

ベクトル(Array)の利用

リストのところで触れたとおり、この演習では配列に当たるものはこのarrayを利用していきます。まず一番簡単な利用方法は以下です。

X = np.array([1, 3, 5, 7])

このarray型のXの値を使いたい場合は、リストと同じようにX[0]などという風に利用できます。

さて、高々数個の配列を作るだけならこれでいいですが、N点ある配列が必要になる場合に、手打ちで入力していては埒が明きません。 C言語では、宣言した後for文で初期化していましたが、ここではNumPyにある関数を利用して、一気に初期化してarray型の配列を作ります。

その方法が以下です。

N = 200
X = np.zeros(N)

これで値が0のN個要素を持つ一次元配列を作ることができ、このXをC言語と同じようにfor文で初期化していくことが出来ます。

例えば、刻み幅がdtが1/Nとなる時間を表すものとし、tにNステップまでの時間を入れておこうとすると以下のようになります。

N = 101
t = np.zeros(N)
dt = 1 / (N-1)
for i in range(N):
    t[i] = i * dt

これで初期化が出来るのですが、これはC言語と同じやり方になります。 しかし、Pythonでは、もっと簡略した方法があります。

N = 101
Tstart = 0
Tend = 1
t = np.linspace(Tstart, Tend, N)

np.linspaceでは、一つ目の引数にスタートの数値、二つ目に終わりの数値、三つ目にその区間を分割する数を入力します。すると、for文を使わなくてもいっぺんに時間の情報を持った配列を準備することが出来ます。

また前に戻りますが、int型で配列を作る場合には、

N = 100
Tstart = 0
Increment = 3
ns1 = np.arange(Tstart, N)
ns2 = np.arange(Tstart, N, Increment)

np.arangeを使えば同様に出来ます。この場合、最後のインクリメントする数を指定しなければ1と自動で指定されます。

小ネタ:リストとnp.array、rangeとnp.arange

リストとnp.arrayはいずれも配列を扱います。ではこれらの間に違いはあるのでしょうか?

list1 = [1, 2, 3]
list2 = np.array([1, 2, 3])

として、変数エクスプローラーで見てみましょう。どちらも[1, 2, 3]という配列ではありますが、それぞれのTypeが、list1では”list”型、list2では”Array of int64″型となっているのが確認できると思います。
まず、Array of int 64型から説明しますと、これは”64 bitのint型の配列”を意味します。すなわち、定義時に全ての構成要素がint型だと、list2の構成要素はすべてint型になるわけです。このときlist2[1]=3.5としてもlist2はint型ですので、自動的にlist2[1]=int(3.5)と同じことが起き、list[1]には3が代入されます。
一方で、list型については、構成要素の型がバラバラでも構いません。例えば、次のような定義も通ります。

list3 = ["計算数学演習", 123, 4.56]

もちろん、list3[1] = ["お腹すいた"]のように、各要素の型が異なる代入も可能です。

どちらも一長一短ありますが、例えば配列を使って辞書的な利用や文字列の解析に用いる際はlist型、数値解析に用いる場合はArray of 〇〇型を使う、というように、適材適所で使いましょう。

次はrangeとnp.rangeについても見てみましょう。

list4 = range(5)
list5 = np.arange(5)

として、変数エクスプローラーで見てみましょう。np.arangeで定義したlist5は上述したArray of 〇〇型なのですが、rangeで定義したlist4はrange型という特殊な型になっています。range型も通常のリストやArrayと同様にlist4[1]などで値を参照できますが、list4[1] = 1のように値を変更することはできません。実はrange型は配列状の変数の集まりではなく、開始点と終点と間隔が定義されているに過ぎないのです。

従って、データとして扱う場合にはnp.rangeを、for文などの単純な繰り返しに使う場合はrangeを使うのが良いでしょう。

課題15

xy平面上の曲線1:y=\exp(x)と 曲線2:y=3+xx>0の範囲で交点を一つもつ。 この交点のx座標を近似的に求めるプログラムを作成せよ。 ただし、正しい解とのずれは0.01以下であることが好ましい。 (ヒント:求め方は自由だが、xを少しずつ増やしていき、曲線1、曲線2の上下が入れ替わる場所を近似的に見つければ良い)


Matplotlib:グラフの描画(次回ちゃんと説明します)

MatplotlibはPythonを利用する大きな理由の一つになる、大変有名な可視化用モジュールになります。NumPyと同じようにモジュールを以下のように読み込みます。

それでは早速単純なsin関数を描いてみたいと思います。

import matplotlib.pyplot as plt
import numpy as np

N = 100
dx = 2.0*np.pi / N
x = np.zeros(N)
y = np.zeros(N)
for i in range(N):
    x[i] = i*dx
    y[i] = np.sin(x[i])

plt.plot(x, y)

(Jupyter Notebookで使う場合は、おまじないとして先頭に%matplotlib inlineを書いてください。)

色々とまだ説明していないコードも書いてありますが、実行してもらうとsin関数が現れてくれると思います。

5行目のnp.piは見ての通りπの数値を与えてくれます。 12行目にあるplt.figure()はそれ自身で、グラフを描画する土台を作ります。実は省略できますが、さまざまな値がデフォルト値のままになるので一般的に利用して各種の値を指定します。figsize=(10,8)と指定することで、グラフを描く前のグラフの大きさを指定します。

最後のplt.show()は実際に目で見るために必要な呼び出しになります。

C言語とgnuplotの組み合わせとの大きな違いは、一つ一つの要素を書き出してグラフ化していたのですが、Pythonでは配列に全てのデータを入れてからグラフを書く形式が主流になります。

次回は、より多くの例題と、いろいろなMatplotlibによる描画を学びます。


確認問題

次回以降の内容では、計算数学演習で解いた以下の問題について、アルゴリズムを理解していることが望ましい。

問題1

100以下の奇数のうち3の倍数でも7の倍数でも無いものを全て求め、 それらを画面に表示し、かつ、それらの和も画面に表示するプログラムを作成せよ。

問題2

リストまたはArrayを利用して、二つの10次元ベクトルu, v の和、差、内積を求めるプログラムを作成せよ。

問題3

次の常微分方程式の初期値問題

    \[\frac{dy}{dt} = -y+\sin(t), y(0) = 1\]

をオイラー法を用いて数値的に解いてグラフ化せよ。 時刻$0 \le t \le 20$の範囲を時間刻み幅0.001(=時間刻み数20001で等分割)、50ステップ毎に値を表示。(Pythonでの描画に関しては次回以降行うので、結果の出力までできれば良い。)

昨年度の計算数学演習のホームページなど参照し、適宜復習してください。 (この演習では、上記問題程度のプログラミングはできるものとして演習を行います。多少の復習はしますが。)


 

 

Department of Mathematical and Life Sciences, Graduate School of Integrated Sciences for Life, Hiroshima University