Pythonで数学を学ぼう: 4月 2018

積分

積分とは、微分の逆を行うことです。
ある式を微分したものを積分すると、元の式に戻ります。（但し微分で省略された部分は除く）
例えば、
5*t**2 + 2*t + 8
を微分すると、
10*t + 2
になりますが、これを積分すると、
5*t**2 + 2*t
になり、元の式（微分で省略された部分は除く）に戻っていることが分かります。

# 微分

from sympy import Symbol, Derivative

t = Symbol('t')

St = 5*t**2 + 2*t + 8

d = Derivative(St, t)

dd = d.doit()

print('微分', dd)

# 積分

from sympy import Integral

i = Integral(dd, t)

print('積分', i.doit())

結果は、以下の様になります。

微分 10*t + 2

積分 5*t**2 + 2*t

微分

Derivativeを使えば微分できます。

from sympy import Symbol, Derivative

t = Symbol('t')

St = 5*t**2 + 2*t + 8

d = Derivative(St, t)

dd = d.doit()

print(dd.subs({t: 2}))

5*t**2 + 2*t + 8を微分すると、

2*5*t + 2となるので、

tに2を代入すると、結果は以下になります。

22

極限

sympyのLimitを使えば極限が求められます。
Sは無限大などの色々な記号の集まりです。

from sympy import Limit, Symbol, S

x = Symbol('x')

l = Limit(1/x, x, S.Infinity)

print(l.doit())

結果は０になります。

上の例の、S.Infinity（無限大）を0にすれば、結果はoo（無限大）となります。

数式をグラフ表示

sympyを使えば、数式をグラフに表示できます。
以下の例は、x**2 + 5x + 4を、xが-10から10の範囲にして出力します。

# グラフを描画

from sympy.plotting import plot

x = Symbol('x')

expr = x**2 + 5*x + 4

plot(expr, (x, -10, 10), title='X', xlabel='x', ylabel='x**2+5x+4')

結果は以下の様になります。

複数の式もグラフ表示できます。

# グラフを描画（２つ）

from sympy.plotting import plot

x = Symbol('x')

expr1 = x**2 + 5*x + 4

expr2 = 6*x + 1

p = plot(expr1, expr2, 1, title='X', xlabel='x', show=False)

p[0].line_color = 'r'

p[1].line_color = 'b'

p.show()

line_colorを設定しているのは、sympyのplotは線の表示色が全て同じになってしまうので、別々の色を設定するために行っています。

結果は以下の様になります。

連立方程式の解

sympyを使えば、連立方程式の解も求められます。

# 連立方程式の解

x = Symbol('x')

y = Symbol('y')

expr1 = 2*x + 2*y - 6

expr2 = 3*x + 2*y - 12

d = solve((expr1, expr2), dict=True)

print(d)

答えは、

[{x: 6, y: -3}]

と表示されます。

解が合ってるか確認するには、subsで上の答えをそれぞれの式に代入してみて0が返却されるか確認します。

# 解が合ってるか検証（0が出力されれば合ってることが分かる）

print(expr1.subs({x:6, y:-3}))

print(expr2.subs({x:6, y:-3}))

解が合ってるので（当然ですが）、下の様に表示されます。

0

0

文字列を数式化する

プログラムで数式を生成するのではなく、文字列から数式を生成するには、sympyのsympifyを使います。

from sympy import sympify

from sympy.core.sympify import SympifyError

expr = (input('Enter a mathematical expression: '))

try:

  expr = sympify(expr)

except SympifyError:

  print('Invalid input')

式の展開

因数分解した後の式を展開するには、Sympyを使います。
例えば、
y=(x+2)(x+3)
の式を展開すると、
y=x**2 + 5*x + 6
になりますが、このように展開式を求めるには、以下の様にします。

>>> import sympy

>>> y = (x + 2) * (x + 3)

>>> y

(x + 2)*(x + 3)

>>> sympy.expand(y)

x**2 + 5*x + 6

因数分解

因数分解をするにはライブラリのSympyを使います。factorで因数分解をします。

>>> import sympy

>>> x = sympy.Symbol('x')

>>> y = sympy.Symbol('y')

>>> y = x**2 + 5*x + 6

>>> y

x**2 + 5*x + 6

>>> sympy.factor(y)

(x + 2)*(x + 3)

相関係数とは、２つの集合がどの程度関係があるかを示す指標となるものです。
例えば、ある月の気温の集合と、同一月のアイスクリームの売り上げの集合があるとして、気温とアイスクリームの売り上げという集合が、どの程度相関しているかを示すのに使われたりします。

相関係数は-1〜1までの範囲で取得でき、-1、または1に近いほど２つの集合が相関していることを示しています。
ただし、相関係数が-1、または1に近いからといって必ず相関しているわけではなく、一応の目安程度のものとなっています。
例えば上の例で言うと、アイスクリームの売り上げが高かった月と同じ月に殺人事件が多かったとしたら、相関係数は多分1に近くなるかもしれませんが、殺人事件とアイスクリームの売り上げが関係しているとはどう考えても思えません。

相関係数は、以下の式で求められます。

２つの集合の要素数：N（２つの集合の要素数は同じである必要があります。）
集合X：x
集合Y：y

NΣxy - ΣxΣy
------------------------------------
√(NΣx²-(Σx)²)(NΣy²-(Σy)²)

def correlation(x, y):

  ''' ２つの集合の分散係数を求める 

      x:集合x

      y:集合y

      集合xと集合yの要素数は同一であること '''

  mulsum = 0

  sum_squared_x = 0

  sum_squared_y = 0

  for xi, yi in zip(x, y):

    # x, yの積和

    mulsum = mulsum + xi * yi

    # xの二乗和

    sum_squared_x = sum_squared_x + xi ** 2

    # yの二乗和

    sum_squared_y = sum_squared_y + yi ** 2

  # xの和

  sum_x = sum(x)

  # yの和

  sum_y = sum(y)

  # xの和の二乗

  squared_sum_x = sum_x**2

  # yの和の二乗

  squared_sum_y = sum_y**2

  N = len(x)

  # 分子

  numerator = N * mulsum - sum_x * sum_y

  # 分母

  denominator = ((N * sum_squared_x - squared_sum_x)

     * (N * sum_squared_y - squared_sum_y)) ** 0.5

  return numerator / denominator

x = [1, 2, 3, 4, 5, 6]

y = [10, 11, 12, 13, 14, 15]

c = correlation(x, y)

print('correlation of x, y is ', c)

中央値

中央値は配列の中央の位置にある要素の値です。
要素数が奇数の場合はぴったり真ん中の位置が求まりますが、偶数の場合は真ん中の要素が２つになるので、真ん中の２つの要素の平均が中央値となります。

def median(arr):

  N = len(arr)

  sorted_arr = sorted(arr)

  if N % 2 == 0:

    # 偶数

    val1 = sorted_arr[int(N / 2) - 1]

    val2 = sorted_arr[int(N / 2)]

    return (val1 + val2) / 2

  else:

    # 奇数

    return sorted_arr[int((N + 1) / 2) - 1]

arr = [3, 6, 1, 10, 8, 7, 9]

med = median(arr)

print('median of array{0} is {1}'.format(arr, med))

median of array[3, 6, 1, 10, 8, 7, 9] is 7

フィボナッチ数列

フィボナッチ数列とは、各要素が前の２つの要素の和となった数列です。
要素数が多くなるにつれて、隣り合う２つの要素の比率が黄金比（1.61803・・・）に近づいていきます。

>>> def fib(n):

...   if n == 1:

...     return [1]

...   elif n == 2:

...     return [1, 1]

...   a = 1

...   b = 1

...   arr = [a, b]

...   for i in range(n):

...     c = a + b

...     arr.append(c)

...     a = b

...     b = c

...   return arr

...

>>> fib(10)

[1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144]

>>> 144/89

1.6179775280898876

最後の２要素の比率が黄金比に近いです。

フィボナッチ数列を１０個しか出していないので少数以下３桁目以降が違いますが、もっと多くの要素数を出すと黄金比（1.61803・・・）に近づいていきます。

投射運動

投射運動とは、例えばボールを投げた時にボールが描く軌跡のことです。
ここでは簡略化のため、水平方向の力にかかる空気抵抗などは省略しています。

ボールが地面に落ちるまでの時間（秒）は、
2 * (初速 * sin(ボールを投げる角度)) / 9.8
で求められます。

from matplotlib import pyplot as plt, axis

import math

def draw_graph(x, y):

  plt.axis(ymax=x[len(x)-1])

  plt.plot(x, y)

def frange(start, end, interval):

  numbers = []

  while start < end:

    numbers.append(start)

    start = start + interval

  return numbers

def draw_trajectory(u, theta):

  theta = math.radians(theta)

  g = 9.8

  # 地面に落ちるまでの時間

  t_flight = 2 * u * math.sin(theta) / g

  intervals = frange(0, t_flight, 0.001)

  x = []

  y = []

  for t in intervals:

    x.append(u * math.cos(theta) * t)

    y.append(u * math.sin(theta) * t - 0.5 * g * t * t)

  draw_graph(x, y)

# 初速25m/sでボールを角度45度で投げた時の軌跡

draw_trajectory(25, 45)

# 初速40m/sでボールを角度45度で投げた時の軌跡

draw_trajectory(40, 45)

# 初速60m/sでボールを角度45度で投げた時の軌跡

draw_trajectory(60, 45)

plt.show()

実行結果は、以下のようになります。

２物質の引力

２つの質量を持つ物質の引力は、下記のようにします。

>>> G = 6.674 * (10 ** -11)

>>> # Gは重力定数

>>> # 距離100m離れた物質、0.5kgと10kgの物質がお互いに引き付け合う引力

>>> force = G * (0.5 * 10.0) / (100.0 ** 2)

>>> force

3.337e-14

>>>

登録: 投稿 (Atom)

Pythonで数学を学ぼう

積分

偏微分

微分

極限