質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Q&A

0回答

2329閲覧

多次元の非線形連立方程式の解き方

aki_ta

総合スコア0

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

0グッド

0クリップ

投稿2020/08/22 07:01

前提・実現したいこと

現在、非線形波動方程式の数値計算をしています。
その際に多次元の非線形連立方程式を解かなければいけません。
そこでpython3のfsolveを使ってコードを書きました。

取り組みの手順として、
まず多次元の線形連立方程式をfsolveで解き、非線形連立方程式を解こうと考えました。
それぞれのコードを書き実行した結果次のようなエラー・警告・問題点が出ました。

・[線形連立方程式]
1.RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
warnings.warn(msg, RuntimeWarning)
2.出力値が違う

・[非線形連立方程式]
TypeError: 'int' object is not callable

・[できれば、、、]
できれば、より早く計算を行いたいと考えています。
※方程式を解けていないため検討をしていません。
しかし、疎行列を扱うためscipy.sparseを用いればよいのではないかと考えています。

それぞれを解決したいです。
下の非線形連立方程式を解ければよいのでfsolve以外でも大丈夫です。
どうぞよろしくお願い致します。
###検討する非線形連立方程式
手書きで見にくいですが解きたいスキームとそれを行列で表した式です。
m:時間空間を分割した際のステップ数
k:x空間を分割した値

イメージ説明
イメージ説明

該当のソースコード

python3

1%matplotlib nbagg 2import numpy as np 3import matplotlib.pyplot as plt 4import matplotlib.animation as animation 5from scipy import optimize 6from scipy.linalg import toeplitz 7from scipy.optimize import fsolve 8 9T=20 10N=2000 11dt=T/N 12L=6 13K=300 14dx=L/K 15r=dt**2/(dx**2) 16R=dt**2/dx 17u=np.zeros((K+1,N,)) 18#二階差分 19D2=toeplitz([-2, 1, *np.zeros(K-1)]) 20 21A=2/r*np.eye(K+1)-D2 22A[0,0]=1+1/R 23A[0,1]=1/R-1 24A[K,K]=1/R+1 25A[K,K-1]=1/R-1 26B=4/r*np.eye(K+1) 27B[0,0]=-1/R 28B[0,1]=-1/R 29B[K,K]=-1/R 30B[K,K-1]=-1/R 31C=-2/r*np.eye(K+1)+D2 32C[0,0]=-1-1/R 33C[0,1]=1-1/R 34C[K,K]=-1-1/R 35C[K,K-1]=1-1/R 36 37#初期値の設定 38for k in range(0,K+1): 39 u[k,0]=np.exp(-(k*dx-L/2)**2) 40 u[k,1]=u[k,0]-dt*2*(k*dx-L/2)**2*np.exp(-(k*dx-L/2)**2) 41 42#----------------------非線形(線形)連立方程式を解く------------------- 43def func(x): 44 x1=u[:,n-1] 45 x2=u[:,n] 46 x3=u[:,n+1] 47#非線形の場合-1/4(x1**3+x3**2*x1+x3*x1**2+x1**3)をつける 48 return np.dot(A,x3)-np.dot(B,x2)-np.dot(C,x1)#-1/4(x1**3+x3**2*x1+x3*x1**2+x1**3) 49#数値計算 50for n in range(1,N-1): 51 x0=u[:,n] 52 u[:,n+1]=fsolve(func, x0) 53#-------------------------------------------------------------------- 54 55# matplotlibで表示 56x = np.linspace(0, L, K+1) # maplotlibで表示するためのx軸の設定 57fig = plt.figure(figsize=(6,4)) 58ax = fig.add_subplot(1,1,1) 59# アニメ更新用の関数 60def update_func(i): 61 # 前のフレームで描画されたグラフを消去 62 ax.clear() 63 ax.plot(x, u[:,i], color='red')#力学的境界条件 64 # 軸の設定 65 ax.set_ylim(-1.2, 1.2) 66 # 軸ラベルの設定 67 ax.set_xlabel('x', fontsize=12) 68 ax.set_ylabel('u', fontsize=12) 69ani = animation.FuncAnimation(fig, update_func, frames=N, interval=1,repeat=False) 70#ani.save("output_riki_3.mp4", writer="ffmpeg", fps=30, bitrate=1000) 71plt.show()

試したこと

1.別の手法を使って線形連立方程式を解き、解があっているかどうかを確かめました。
そのため、上のコードの「非線形(線形)連立方程式をとく」の中身がエラー・警告の原因だと考えています。
下記に線形連立方程式を解いたコードを書きます。
2.インターネット上では3次元の非線形連立方程式をとくサイトはありましたが多次元のサイトはなく、見よう見まねで先ほど書いたコードを書きました。

python3

1%matplotlib nbagg 2import numpy as np 3import matplotlib.pyplot as plt 4import matplotlib.animation as animation 5from scipy import optimize 6from scipy.linalg import toeplitz 7import scipy.sparse 8import scipy.sparse.linalg as spla 9 10T=20 11N=2000 12dt=T/N 13L=6 14K=300 15dx=L/K 16r=dt**2/(dx**2) 17R=dt**2/dx 18u=np.zeros((K+1,N,)) 19 20A=2/r*np.eye(K+1)-D2 21A[0,0]=1+1/R 22A[0,1]=1/R-1 23A[K,K]=1/R+1 24A[K,K-1]=1/R-1 25B=4/r*np.eye(K+1) 26B[0,0]=-1/R 27B[0,1]=-1/R 28B[K,K]=-1/R 29B[K,K-1]=-1/R 30C=-2/r*np.eye(K+1)+D2 31C[0,0]=-1-1/R 32C[0,1]=1-1/R 33C[K,K]=-1-1/R 34C[K,K-1]=1-1/R 35 36#初期値の設定 37for k in range(0,K+1): 38 u[k,0]=np.exp(-(k*dx-L/2)**2) 39 u[k,1]=u[k,0]-dt*2*(k*dx-L/2)*np.exp(-(k*dx-L/2)**2) 40 41#数値計算 42for n in range(1,N-1): 43 b=np.dot(B,u[:,n])+np.dot(C,u[:,n-1]) 44 a_csr = scipy.sparse.csr_matrix(A) 45 u[:,n+1] = spla.dsolve.spsolve(a_csr, b) 46 47# matplotlibで表示 48x = np.linspace(0, L, K+1) # maplotlibで表示するためのx軸の設定 49fig = plt.figure(figsize=(6,4)) 50ax = fig.add_subplot(1,1,1) 51# アニメ更新用の関数 52def update_func(i): 53 # 前のフレームで描画されたグラフを消去 54 ax.clear() 55 ax.plot(x, u[:,i], color='red') 56 #軸の設定 57 ax.set_ylim(-1.2, 1.2) 58 # 軸ラベルの設定 59 ax.set_xlabel('x', fontsize=12) 60 ax.set_ylabel('u', fontsize=12) 61ani = animation.FuncAnimation(fig, update_func, frames=N, interval=1,repeat=False) 62#ani.save("output_riki_3.mp4", writer="ffmpeg", fps=30, bitrate=1000) 63plt.show()

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

まだ回答がついていません

会員登録して回答してみよう

アカウントをお持ちの方は

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問