min117の日記

初期desireもち。趣味Mac,メインFedora,仕事xp。

Python数学 y=a**x と y=b**x の交わるトコ(積分の面積部分)に色を付ける

Pythonで、このグラフの斜線に色つけられるかなー、とか。

youtu.be

動画のなかで出てきた式は違うよ、ってChatGPTさんは言うのだが…

ま、いいや。コードしてみる。

y = a**x と y = -x**2 + b で囲まれた部分の面積

できた。

もうね、ChatGPTあれば何でもできるね。

 

$ cat myKINRImath19d.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.optimize import fsolve

#import pandas as pd
import pandas_datareader as web
#import matplotlib.pyplot as plt
# %matplotlib inline
import matplotlib.font_manager as fm
# plt.rcParams['font.family']='Droid Sans Japanese'
import japanize_matplotlib

from matplotlib import rcParams
rcParams['figure.figsize']=15,10
rcParams['font.size']=15
import sys
import os
import math

#import cv2
#from PIL import Image, ImageFont, ImageDraw
from PIL import Image
from PIL import ImageDraw
from PIL import ImageFont
import numpy as np

import subprocess

from sympy import Symbol, Derivative, sympify, pprint, pretty
from sympy.core.sympify import SympifyError
# ChatGPT from sympy import sympify, Symbol

import datetime
from datetime import datetime as dt
# import sys
# import os
# import numpy as np
# import matplotlib.pyplot as plt

import subprocess
# cdしないとcronで/rootに溜まっちゃう!
# スクリプトディレクトリを取得
script_directory = os.path.dirname(os.path.abspath(__file__))

# スクリプトディレクトリにカレントディレクトリを変更
os.chdir(script_directory)
#os.chdir('/media/WD30EZRX/PT3/')

# 現在の日時を取得
#execution_time = datetime.now().strftime("%Y-%m-%d %H:%M:%S")
myTIMEE = datetime.datetime.now().strftime('%Y/%m/%d %H:%M:%S')
myARG0 = os.path.basename(__file__)

# myARG1 = sys.argv[1]


mySOURCE = ('''積分の面積公式5選!一瞬で面積を求める裏ワザを教えます\n https://www.youtube.com/watch?v=gPgGd22Nq0o''')

#mySOURCE = ('''
#ネイピア数 自然対数の底e とは\n
#    https://www.youtube.com/watch?v=1M7FF1nd25I\n
#に出てきたグラフを描く''')
print(mySOURCE)
print('-------------------------------')

 

 

# 定数
a = 2  # 例: a^x の a の値
b = 4  # 例: b^x の b の値
# alpha = 0  # 例: x = α
#alpha = -2  # 例: x = α

print('a を浮動小数点数に変更')
a = 2.0  # 'a' を浮動小数点数に変更
beta = 0  # 例: x = β

print('関数定義')
def f_a(x):
    return a ** x

def f_b(x):
    #return b ** x
    return -x**2 + b

print('交点を探すための関数')
def intersection(x):
    return f_a(x) - f_b(x)

print('交点を見つける(2つの交点を fsolve で探す)<- new')
alpha, beta = fsolve(intersection, [-2, 1])
print(f"交点: alpha={alpha:.4f}, beta={beta:.4f}")

print('面積計算の関数')
def integrand(x):
    # return np.abs(f_a(x) - f_b(x))
    return f_a(x) - f_b(x)  # f_a(x) が常に大きい場合← f_b(x)が下に凸

area, _ = quad(integrand, alpha, beta)
print(f"面積: {area:.4f}")

print('xの範囲')
# x = np.linspace(alpha, beta, 500)
# x = np.linspace(-4, 4, 1000)
x = np.linspace(alpha - 2, beta + 2, 1000)


print('グラフの描画')
#OKK plt.figure(figsize=(8, 6))
#` plt.plot(x, f_a(x), label=r'$y=a^x$', color='blue')
#` plt.plot(x, f_b(x), label=r'$y=b^x$', color='red')
plt.plot(x, f_a(x), label=r'$y=a^x$', color='blue', linewidth=2)
plt.plot(x, f_b(x), label=r'$y=-x^2 + b$', color='red', linewidth=2)

#OKK # 塗りつぶし
#OKK plt.fill_between(x, f_a(x), f_b(x), where=(f_a(x) > f_b(x)), color='lightblue', alpha=0.5)
#OKK plt.fill_between(x, f_b(x), f_a(x), where=(f_b(x) > f_a(x)), color='lightcoral', alpha=0.5)

# # 塗りつぶし
# plt.fill_between(x, f_a(x), f_b(x), color='lightblue', alpha=0.5)

#@ print('塗りつぶしぃ')
#@ plt.fill_between(x, f_a(x), f_b(x), where=(f_a(x) > f_b(x)), color='lightblue', alpha=0.5)

#OKK print('塗りつぶし: 交わる点と x=0 の間で塗りつぶし')
#OKK x_fill = np.linspace(alpha, beta, 500)
#OKK plt.fill_between(x_fill, f_a(x_fill), f_b(x_fill), color='lightgreen', alpha=0.6)

print('交点間の塗りつぶし')
x_fill = np.linspace(alpha, beta, 500)
plt.fill_between(x_fill, f_a(x_fill), f_b(x_fill), color='lightgreen', alpha=0.6)


print('軸と凡例')
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)

print('x軸とy軸の範囲を調整 <-new')
plt.xlim(alpha - 2, beta + 2)
plt.ylim(min(f_b(x)), max(f_a(x)) + 1)

#OKK print('x軸とy軸の範囲を広く設定')
#OKK plt.xlim(-4, 4)
#OKK plt.ylim(-10, 10)

# plt.title(r'Area between $y=a^x$ and $y=b^x$')
plt.title(r'Area between $y=a^x$ and $y=-x^2 + b$ (between intersections)')

# plt.xlabel('x')
# plt.ylabel('y')
plt.xlabel('x', fontsize=14)
plt.ylabel('y', fontsize=14)
plt.legend()

print('グラフ表示')
plt.grid(True)
plt.legend()
#plt.show()

# コメントをx軸の下に追加
# plt.figtext(0.5, -0.1, f"Script: {myARG0} | Executed on: {execution_time}", ha="center", fontsize=10)
# plt.figtext(0.5, -0.1, f"Script: {myARG0} | Executed on: {myTIMEE}", ha="center", fontsize=10)
kawacom = ''
#kawacom = f'{myARG1}のグラフと接線\n{datetime.datetime.now().strftime("%Y/%m/%d %H:%M:%S")} \n({myARG0})'
#kawacom = f'{datetime.datetime.now().strftime("%Y/%m/%d %H:%M:%S")} \n({myARG0})'
kawacom = f'{datetime.datetime.now().strftime("%Y/%m/%d %H:%M:%S")} ({myARG0}) \n' + mySOURCE
plt.figtext(0.5, 0.01, kawacom, ha="center", fontsize=10)

#kawacom =  datetime.datetime.now().strftime('%Y/%m/%d %H:%M:%S') + ' \n(' + myARG0 + ')'
#@ kawacom = mySOURCE + '\n' + datetime.datetime.now().strftime('%Y/%m/%d %H:%M:%S') + ' \n(' + myARG0 + ')'

## plt.xlabel(kawacom)
# コメントをx軸の欄外に追加
# plt.figtext(0.5, -0.2, kawacom, ha="center", fontsize=10)
plt.figtext(0.5, 0.01, kawacom, ha="center", fontsize=10)

# 数式を吹き出しコメントとして表示
# plt.annotate(rf"${equation_input}$", xy=(0.5, 0.5), xytext=(0.5, 0.75),
#     textcoords='axes fraction', ha='center', va='center',
#     bbox=dict(boxstyle='round,pad=0.5', fc='yellow', alpha=0.5),
#     arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0.5'))


#math12 # グラフをプロット
#math12 plt.figure(figsize=(8, 6))
#math12 plt.plot(x, y, label=r"$x^5 - 30x^3 + 50x$")
#math12 plt.xlabel('x')
#math12 plt.ylabel('y')
#math12 plt.title('Plot of $x^5 - 30x^3 + 50x$')
#math12 plt.grid(True)
#math12 plt.legend()

myPNG = 'myKINRImath19d.png'

# PNG形式で保存
# plt.savefig('myKINRImath12.png')
# plt.savefig(myPNG)
plt.savefig(myPNG, bbox_inches='tight')
#plt.show()
print('-------------------------------')
print(myPNG + ' を保存しました')
#print('(yの座標軸を画面中央に表示した版)')
print('-------------------------------')
print(mySOURCE)
print('-------------------------------')
print('''
積分の面積公式5選!一瞬で面積を求める裏ワザを教えます
\n
https://www.youtube.com/watch?v=gPgGd22Nq0o
''')
print('-------------------------------')
print('http://192.168.3.11:5000/PT3/doLSglob.php?CHDIR=PT3&GREPWORD=myKINRI&MKSORT=goMKSORT')