Python 学習:複雑?な方程式を解く

2015-12-18 :  PCクリニック

本文の前に、
-・・・ -・-
現時点での blogramのランクインカテゴリは、
6、3、0、0、1、 0、0、0、0、0(41)で、換算ポイント 86pt 。
「Perl」「化学業界」「硝子業界」「e-radio」「グルコサミン」、
「Python」「FM COCOLO」「Firefox」bg値変動のみ。
「C言語」「FM青森」全く変化無し。
・-・ - -・

さて、本文。

先日(2015-12-12)の記事:「積分関数を求めるには(不定積分)
では、“Python”は“Mathematica”に敵わなかった?

そこで、(<紙>にとっては)複雑な方程式を解く問題で比較してみた。

対象となるのは、
某学会誌に載っていた方程式( から<紙>流に書き直したもの )。
(* 関数4つ + 1つ *)
Ge[x_] := 10^(A12*(A21*(1-x)/(A12*x+A21*(1-x)))^2);
Gw[x_] := 10^(A21*(A12*x/(A12*x+A21*(1-x)))^2);
Je[x_, y_] := Pe*(Ge[x]*x-Ge[y]*y*(Exp[Ve*(P2-P1)/(R*T)]));
Jw[x_, y_] := Pw*(Gw[x]*(1-x)-Gw[y]*(1-y)*(Exp[Vw*(P2-P1)/(R*T)]));
F[y_] := Je[x, y]*(1-y)-Jw[x, y]*y;

(* 変数に値設定 *)
A12 = 0.7292; A21 = 0.4104;
Pe = 0.01; Pw = Pe*10; Ve = 5.85*10^(-5); Vw = 1.8*10^(-5);
P1 = 30.0*10^6; P2 = 0.1*10^6; R = 8.31451; T = 300.0;

(* 解きたい方程式は、変数 x に対して F[y]==0 となる y の値 *)
x = 0.3; FindRoot[ F[y] == 0 , {y, 0, 1}]
実行すると、
  { y -> 0.190151 }
と出る。


そこで、
これを“Python”で書き直した。
使ったのは、先月(2015-11-27)の記事:
Python 学習:scipy.optimize.fsolve
に従った。
# -*- coding: utf-8 -*-
import numpy as np
import scipy.optimize

### 関数4つ+1つ
def Ge(x): return 10**( A12*( A21*(1-x) / ( A12*x + A21*(1-x) ) )**2)
def Gw(x): return 10**( A21*( A12*x / ( A12*x + A21*(1-x) ) )**2)
def Je(x,y): return Pe*( Ge(x)*x - Ge(y)*y*( np.exp(Ve*(P2-P1)/(R*T))))
def Jw(x,y): return Pw*( Gw(x)*(1-x) - Gw(y)*(1-y)*( np.exp(Vw*(P2-P1)/(R*T))))
def F(y): return Je(x,y)*(1-y) - Jw(x,y)*y

### 変数に値設定
A12=0.7292; A21=0.4104
Pe=0.01; Pw=Pe*10; Ve=5.85*10**(-5); Vw=1.8*10**(-5)
P1=30.0*10**6; P2=0.1*10**6; R=8.31451; T=300.0

### 解きたい方程式は、変数 x に対して F(y)==0 となる y の値
x=0.3; y=scipy.optimize.fsolve( F, [0] )
print y
実行すると、
  [ 0.19015104 ]
と出た。


この位だと、それほど複雑では無いのかナ?

安心した。


本日はここまで。


見ていただいた序でとは厚かましい限りですが、
お帰りに投票して頂けるとなお嬉しいです。 ⇒ blogram投票ボタン


151212
関連記事
スポンサーサイト

コメントの投稿

管理者にだけ表示を許可する

おきてがみ/blogram
blogram投票ボタン



おきてがみ

最新記事
カレンダー
07 | 2017/08 | 09
- - 1 2 3 4 5
6 7 8 9 10 11 12
13 14 15 16 17 18 19
20 21 22 23 24 25 26
27 28 29 30 31 - -
月別アーカイブ
カテゴリ
最新コメント
検索フォーム
リンク
プロフィール

<紙>

Author:<紙>
ようこそ。
「パソコンヲタクの雑記帳」
もろもろなことを綴っています。
パソコン ヲタクってねくら?
画像は kami でなく kani です。

カウンター(fc2、i2i) /Google Analytics


i2i(from 2010-08-24)
Total =
Today  =  
Yesterday=
アンチエイジング

Google Analytics
ブックマーク