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


第１・２回（第一週）

はじめに

計算数理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

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

ここに3とありますが、pythonには現在2系と3系と呼ばれる異なるバージョンが広まっています。しかし将来的には3系のみになる予定なので、ここにあるようにpython3と実行するようにしてください。

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

#include <stdio.h>

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

    return 0;
}

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

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

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

  • ;（セミコロン）がpythonでは書いてはいけない。
  • main関数が存在しない。

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

Jupyterを使う

さらに大きな武器を使ってプログラムの学習を進めていきたいと思います。
これまではエディタでプログラムを書き、terminalでコンパイルして実行してファイルにデータを保存して、gnuplotでグラフを書きました。

Pythonではこれらの処理を一つのブラウザ上の画面で行います。
数式処理演習で用いたmathematicaととてもよく似ていると思います。
しかし、言語としてはPythonは言語の書き方としてC言語にかなり近いので、mathematicaとはかなり違うので注意してください。

それでは、Jupyterを実行しましょう。まずターミナルを開いて、

jupyter notebook

と実行してください。しばらくするとデフォルトブラウザで次の様な画面が現れるはずです。

ここにいろいろなファイル・ディレクトリが出ていますが、ターミナルを開いてすぐにlsを実行したときに見られるものが見られているはずです。

次に、右上にあるNewを選択して、

Python3というものをクリックしてください。

すると最終的に、次の様な画面が出ます。

これでJupyterのノートブックというファイルが作られます。
このファイル名は、Untitledになっているので、その部分をクリックすると以下のような画面になり、ファイル名を変更できます。

今回はLec1としてRenameをクリックしてファイル名を変更します。

これで利用環境は整いました。
Jupyterでは灰色の部分にプログラムを書き込みます。この部分をCellと呼びます。その部分にさきほどのHello, World!!を書き出すプログラムを書き込み、そのCellを選択した状態（青い枠で囲われている状態）でRunというボタンを押すとそのすぐ下に出力されます。

データ型・変数・初期化

それでは、Pythonの内容に進んでいきますが、C言語との大きな違いが変数にも表れます。

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

C言語ではint型やfloat型の変数の宣言について学びました。
Pythonにも同様に型は存在し、データ型と呼ばれていますが、変数の宣言という

float x;

のような形で変数は宣言できません。
その代わりに宣言と初期化を以下のように行います。

x = 0.0

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

x = 0

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

a = 1
x = 2.1

print(a, x)

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

リスト

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

X = [1, 3, 5, 7]

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

print(X[1])

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

X[1] = 0

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

制御文：if文

次にif文について学びます。これはC言語と少し異なるのと、pythonを学ぶ上で重要な「インデント」について学びます。
C言語ではif文で条件式は（）に囲まれ、真の場合に実行されるプログラムは、{}に囲まれていました。
Pythonでは囲まれず、条件式は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かつｂが1
or || if a==0 or b==1: aが0またはｂが1
not ! if not a==1: aが0ではない

課題１

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

Pythonの入力について

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

s = input()

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

s = input()

# int型
s = int(s)

# float型
s = float(s)

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

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

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

課題２

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

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

課題３

２つの実数 s, t を入力すると, s ＞ t, s = t, s ＜ t か判定して表示するプログラムを作成せよ．（If-else文のパターン３を参考）

課題４

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

課題５

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

課題6

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

制御文：for文・while文

次に重要なfor文についてです。これはC言語とかなり異なります。
C言語のfor文では初期化式・条件式・反復式の三つがありました。
pythonでは、0から9まで出力するには以下のようになります。

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

とするだけで繰り返しが実行されます。
このrangeというものはいくつか使い方があり、このように一つだけ数字を入れると、0からその数まで数を繰り返していきます。

N回繰り返して計算したい場合には、以下のようになります。

N = 100

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

rangeには次の様な使い方があります。

N = 100

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

これを実行すると1, 4, 7, …, 97という数字が出力されるようになります。
何をしているかというと，rangeの引数にNの前に1，後に3が加わることで，まず初めの1がfor文のスタート地点，最後の3が増加する量を指定します。
C言語で考えると

int i;
int N = 100;

for (i=0; 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文はインデントで制御文の範囲を指定する以外は、C言語に近いです。

例えば、

N = 100
i = 0

while i < 20:
    i = i + 1 
    print(i)

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

これだけでは機能が少なくPythonを学ぶありがたみが少ないです。
Pythonを利用するありがたみは、その豊富なライブラリにあります。
（もちろん他にもありますが。。。）
本講義で主に用いる、NumPyとMatplotlibというライブラリに使用方法を学び、実際に問題を解いてpythonに慣れていきたいと思います。

それでは，例題に取り組んでいきましょう。

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

wa = 0
seki = 1

for i in range(5):
    wa = wa + i
    seki = seki * i;

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

課題１

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

課題２

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

課題３

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

課題４

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

課題５

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

以下，難しめの問題

課題６

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

課題７

7のn乗をAnとあらわすとする．下三桁が(943)となる最小のAnのnを求めよ．
(追記：nが存在しない場合は，そのように出力せよ．）
（intで扱える整数は，-2147483648から2147483647であることに注意．変数の値が2147483647を超えないように工夫しつつ，下三桁の値を調べていくこと．）

課題8

整数Mのn乗をMnとあらわすとする．三桁以下の整数Mの入力を促し下二桁が(16)となる最小のMnのnを求めよ.
また，そのようなnが存在しない場合は，「そのようなnは存在しません」と出力すること．

課題9

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

NumPy

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

import numpy as np

C言語で#includeに対応するものがimoprtになります。そのまま後ろにas npとありますが，詳細は触れませんがプログラムの中ではこのnpというものを多用します。
importの後ろに書くものは，一般的にライブラリではなくモジュールと呼ばれます。特にNumPyはNumericalのNumから名前の由来のある数値計算用のモジュールで，これがあればかなりのことは出来ます。

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となる時間を表すものとし，XにNステップまでの時間を入れておこうとすると以下のようになります。

N = 100
dt = 1 / N

for i in range(N):
    x[i] = i * dt

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

N = 100
Tstart = 0
Tend = 2

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

 

N = 100
Tstart = 0
Tend = 2

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

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

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

N = 100
Tstart = 0
Increment = 3

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

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

Matplotlib

Matplotlibはpythonを利用する大きな理由の一つになる、大変有名な可視化用モジュールになります。NumPyと同じようにモジュールを以下のように読み込みます。ただし，Jupyterで使う場合のおまじないがありますが，セットで書いておいてください。

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

%matplotlib inline
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.figure(figsize=(10,8))
plt.plot(x, y)
plt.show()

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

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

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

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

来週は，より多くの例題と，いろいろなMatplotlibによる描画を学びます。

関数

関数を利用者が自ら定義して利用することは，どのプログラミング言語でも必要な基礎的な技術です。早速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言語と違い先に関数を宣言しておくプロトタイプ宣言ができません。基本的に宣言という概念がないので当然ではあるのですが，常に利用するよりも前に関数を定義しておく必要があるので注意してください。

課題

最終的に以下の問題を解けるようになってから，演習を進めていけるようにしたいと思います。

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

問題2
配列変数を利用して、二つの10次元ベクトルu, v の内積を求めるプログラムを作成せよ。

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

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


をオイラー法を用いて数値的に解いてグラフ化せよ。
(時刻$0 \le t \le 20$の範囲を20000等分割、50ステップ毎に値を表示)

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