fortran66のブログ

fortran について書きます。

【メモ帳】積分2

C6.1 \int_0^{\pi\over2} \cos^{-1}\left( {1\over1+2\cos x} \right)  dx=\zeta(2)

I=\int_0^{\pi\over2} \arccos \left( {1\over1+2\cos x} \right)  dx

ここで  \cos2x=2\cos^2x-1=2u^2-1, u=\cos x より、
\arccos(2u^2-1)=\arccos(\cos2x) = 2x = 2\arccos(u) なので、
\alpha=2\theta^2-1 とすると、\theta=\sqrt{\alpha\over 2} すなわち \arccos(\alpha)=2\arccos\sqrt{1+\alpha\over2} となる。

これより

\begin{align}
\alpha&=\arccos\left({\cos x\over1+2\cos x}\right)\\
          &=2\arccos\left(\sqrt{1+\cos x\over1+2\cos x}\right)\\
&=2\arctan\left(\sqrt{\cos x\over1+\cos x}\right)\\
\end{align}
が得られる。ここで最後の式はピタゴラスの定理三角関数の定義による。

 
\begin{align}
I&=2\int_0^{\pi\over2}\arctan\left( \sqrt{\cos x\over 1+\cos x}\right) dx\\
&=4\int_0^{\pi\over4}\arctan\left( \sqrt{\cos2y\over1+\cos2y} \right)dy, \,\,(x=2y)\\
&=4\int_0^{\pi\over4}\arctan\left( \sqrt{1-2\sin^2y\over2\cos^2y} \right)dy, \,\,(倍角の公式)\\
\end{align}

ここで {1\over b}\arctan b =\int_0^1{1\over1+b^2t^2}dt より、
 
\begin{align}
I&=4\int_0^{\pi\over4}\int_0^1\sqrt{1-2\sin^2y\over2\cos^2y} {1\over 1+\left({1-2\sin^2y\over2\cos^2y}\right)t^2} dt dy\\
&=4\int_0^{\pi\over4}\int_0^1{\sqrt2\sqrt{1-2\sin^2y}\,\cos y\over (2+t^2)-2(1+t^2)\sin^2y} dt dy\\
&=4\int_0^{\pi\over2}\int_0^1{\sqrt{1-\sin^2w}\,\cos w\over (2+t^2)-(1+t^2)\sin^2w} dt dw, \,\,(\sqrt2\sin y=\sin w )\\
&=4\int_0^{\pi\over2}\int_0^1{\cos w\cos w\over (2+t^2)-(1+t^2)(1-\cos^2w)} dt dw\\
&=4\int_0^{\pi\over2}\int_0^1{\cos^2w\over 1+(1+t^2)\cos^2w} dt dw\\
&=4\int_0^\infty\int_0^1{{1\over 1+s^2}\over 1+(1+t^2){1\over1+s^2}}{1\over1+s^2} dt ds, \,\,(s=\tan w )\\
&=4\int_0^\infty\int_0^1{1\over (1+s^2)(2+t^2+s^2)} dt ds\\
&=4\int_0^\infty\int_0^1{1\over(1+t^2)(1+s^2)}-{1\over(1+t^2)(2+t^2+s^2)}dtds, \,\, (部分分数に展開)\\
&=4\left({\pi^2\over8}-{\pi\over2}\int_0^1{1\over(1+t^2)\sqrt{2+t^2}} \right) , \,\, ( s について積分ののち t について積分)\\
&=4\left({\pi^2\over8}-{\pi\over2}{\pi\over6} \right), \,\, (昨日の結果から)\\
&={\pi^2\over6}\\
\end{align}

数値計算 台形公式+Richardson 補外

実行結果
   1.64492430955605        1.64490909118690        1.64486343607944
   1.64493406684823
続行するには何かキーを押してください . . .

収束が悪い。

プログラム
module m_mod
    implicit none
    integer, parameter :: kd = kind(0.0d0)
    real(kd), parameter :: pi = 4 * atan(1.0_kd)
contains
    pure real(kd) function f(x)
        real(kd), intent(in) :: x
        f = acos(1.0_kd / (1.0_kd + 2.0_kd * cos(x)))
    end function f
    
    pure real(kd) function trapez(func, x0, x1, n) result(s)
        procedure(f) :: func
        real(kd), intent(in) :: x0, x1
        integer , intent(in) :: n
        real(kd) :: h, y(0:n)
        integer :: i
        h = (x1 - x0) / n
        forall(i = 0:n) y(i) = func(x0 + i * h)
        s = h * ( y(0) / 2 + sum(y(1:n-1)) + y(n) / 2 )
    end function trapez
end module m_mod

program test
    use m_mod
    implicit none
    integer, parameter :: n = 2**9
    real(kd) :: s0, s1, x0, x1
    integer :: i
    x0 = 0.0_kd 
    x1 = pi / 2.0_kd - epsilon(x1) 
    s0 = trapez(f, x0, x1, n * 2)
    s1 = trapez(f, x0, x1, n)
    print *, (4 * s0 - s1) / 3, s0, s1, pi**2 / 6.0_kd
end program test  

【メモ帳】積分


\begin{align}
I&=\int_0^1{dt\over(t^2+1)\sqrt{t^2+2}}\\
&(t=\sqrt2\tan x, dt=\sqrt2(1+\tan^2x)dx, 1=\sqrt2\tan a)\\
&=\int_0^a{\sqrt2(1+\tan^2x)\over(2\tan^2x+1)\sqrt{2\tan^2x+2}}dx\\
&=\int_0^a{{\sqrt2\over\cos^2x}\over\left({2\over\cos^2x}-1\right){\sqrt2\over\cos x} }dx\\
&=\int_0^a{\cos x\over2-\cos^2x}dx\\
&=\int_0^a{\cos x\over1+\sin^2x}dx\\
\end{align}

ここで \tan y = \sin x, [ y=\arctan(\sin x) ]とおくと、

\begin{align}
(1+\tan^2y)dy&=\cos x dx\\
(1+\sin^2x)dy &=\cos x dx\\
dy&={\cos x\over 1+\sin^2x}dx
\end{align}
となるので、求めたい積分I=\int_0^b dy=b となる。ただし  \tan b = \sin a なので

\begin{align}
1&=\sqrt2\tan a\\
{1\over\sqrt2}&={\sin a \over\cos a}\\
{1\over2}&={\sin^2a\over(1-\sin^2a)}\\
1-\sin^2a&=2\sin^2a\\
3\sin^2a&=1\\
\sin a&={1\over\sqrt3} \,\,(>0)\\
\end{align}
よって、
I=b=\arctan(\sin a)=\arctan{1\over\sqrt3}={\pi\over6}
となる。

数値積分 台形公式+Richardson 補外

実行結果
  0.5235988      0.5235969      0.5235909     pi/6=  0.5235988
続行するには何かキーを押してください . . .

プログラム

module m_mod
    implicit none
    real, parameter :: pi = 4 * atan(1.0)
contains
    pure real function f(x)
        real, intent(in) :: x
        f = 1.0 /((x*x+1.0)*sqrt(x*x+2.0))
    end function f

    pure real function trapez(func, x0, x1, n) result(s)
        procedure(f) :: func
        real   , intent(in) :: x0, x1
        integer, intent(in) :: n
        real :: h, y(0:n)
        integer :: i
        h = (x1 - x0) / n
        forall(i = 0:n) y(i) = func(x0 + i * h)
        s = h * ( y(0) / 2 + sum(y(1:n-1)) + y(n) / 2 )
    end function trapez
end module m_mod

program test
    use m_mod
    implicit none
    integer, parameter :: n = 2**6
    real :: s0, s1, x0, x1
    x0 = 0.0 
    x1 = 1.0  
    s0 = trapez(f, x0, x1, n * 2)
    s1 = trapez(f, x0, x1, n    )
    print *, (4 * s0 - s1) / 3, s0, s1, 'pi/6=', pi / 6.0
end program test  

(微修正H30.6.5)

Modern Fortran Explained (Numerical Mathematics and Scientific Computation)

Modern Fortran Explained (Numerical Mathematics and Scientific Computation)

  • 作者: Michael Metcalf,John Reid,Malcolm Cohen
  • 出版社/メーカー: Oxford University Press, U.S.A.
  • 発売日: 2011/05/08
  • メディア: ペーパーバック
  • この商品を含むブログを見る

【寝言】6月4日 支那天安門事件

www.sankei.com
www.sankei.com


政権の正統性を持たない中共や韓国・北朝鮮の如きみじめな未開人は、実に哀れである。18世紀啓蒙思想の水準に達していない。彼らの愚行のニュースを見るたびに苦笑いが浮かぶ。

中世以前の奴隷制社会の金氏朝鮮、選挙の真似事をして見せても、所詮は土人のヒステリー祭りになる韓国や、それらの一匹半のシャム双生児の生みの親であるの母犬の中共。全く救いのない無間地獄を見るようである。

これら土人国に共通するのは、すべてが独裁政治に従属していることで、歴史も文化も芸術も全て今現在の政治に従属する道具にすぎず、客観性や事実性・論理的整合性が欠落しているところにある。

それにつけても、(第二)天安門事件の死亡者は0人だったのだからw堂々と事実を出せばよいw 真実を知りたければ、天安門広場の壁新聞を読めといって、集まった連中を、また戦車で踏み潰せばいいだろう。

かつて支那全土を支配したモンゴル人や色目人の人口比は約1~2%、清朝満州族もその程度、我が国が支那に正義と秩序をもたらしたときはさらに少ない人数だった。

今日の中国共産党の党員数も人口のほぼ1~2%なので、まぁ何も変わらぬ進歩のない土人国である。はやくマルクス主義を捨てて、男は辮髪・女は纏足・学生は四書五経の丸暗記に戻って欲しいものである。

あと糖質制限ダイエットならぬ、僻穀して仙人を目指してください。

狩野直喜の「中国哲学史」(岩波書店なので題名が本来の支那哲学史から改竄されている)によると、すでに漢初の賈誼が胎教を唱えているそうで、淫靡な鄭声でも聞かせて匹夫・凡夫のラッパー目指せよと思います。
www.sankei.com

中国政府は23日までに、ヒップホップ文化は低俗だとして、ラップ歌手らをテレビやラジオ番組に出演させない方針を示した。中国メディアが伝えた。人気曲の薬物使用に関する歌詞が問題視されたことが直接のきっかけだが、体制批判に結び付きやすいヒップホップ文化が大衆に浸透することを警戒したとみられる。

 中国メディアによると、メディアを管理する国家新聞出版ラジオ映画テレビ総局が最近「入れ墨のある芸能人、(ラップを含む)ヒップホップ文化、不健全な文化は番組で扱わない」よう関係機関に求めた。

 中国では昨年末ごろから、人気ラップ歌手の男性の曲に「白い粉末を並べておく」「恥知らずのくそ女」など薬物使用や女性蔑視の表現があるとしてネット上で攻撃対象となった。

 今年に入り、共産党機関紙、人民日報や国営メディアも「青少年に悪影響を与えている」と批判。男性は謝罪し、作品のネット配信も停止された。関係者によると、当局の意向を受け、ラップ歌手の出演を見合わせるテレビ局も出てきている。(共同)

阿片中毒国で文身した匹夫凡夫禁止令w


中国哲学史 (1953年)

中国哲学史 (1953年)

しかし、糞みたいな新書本や文庫本が1,000円近くするというのに、ケース入りの狩野直喜中国哲学史が700円しないとか笑える。

【朗報】Intel Fortran v.18 Update 3 で F1 HELP 復活・・・しかし

Intel Compiler v.18.0.3 出る!

Intel Fortran v.18 update 3 が出ました。

リリースノートによると、今回の機能拡張は以下の通り。

Changes in Update 3 (Intel® Fortran Compiler 18.0.3)

  • Changes to mitigate speculative execution side-channel issues
  • Restored context sensitive help in VS
  • /Q[a]xIcelake-server and /Q[a]xIcelake-client options added to support Ice Lake microarchitecture
  • /Q[a]xCannonlake option added to support Cannon Lake microarchitecture
  • Corrections to reported problems

Intel® Visual Fortran Compiler 18.0 for Windows* Release Notes for Intel® Parallel Studio XE 2018 | Intel® Software

第二項目に Visual Studio での F1 HELP の復活があります。しかしながら、日本語環境下では、Version 17 の Help ファイルの内容が表示されてしまいます。

Option-Environment-International Settings=>Language -> English でメニュー等を英語モードに変更後 VS 再起動で v.18 用の Help ファイルが表示されます。

ただし v.18 用の Help ファイルは手動でセットする必要があります。詳しくは以下のページ参照。
Download Documentation: Intel® Compiler (Current and Previous) | Intel® Software

Fortran ハンドブック

Fortran ハンドブック

Fortran90/95プログラミング

Fortran90/95プログラミング

建築構造設計・解析入門 Fortran 解析プログラム付

建築構造設計・解析入門 Fortran 解析プログラム付

第3版 有限要素法による流れのシミュレーション OpenMPに基づくFortranソースコード付

第3版 有限要素法による流れのシミュレーション OpenMPに基づくFortranソースコード付

【世俗ニュース】言論への忌まわしくおそろしいテロ!

NHK は自局の言論を守るためにも、もっと大々的に報道すべき!

NHK近くで制作会社員切りつけた疑い 韓国籍の男逮捕

www.asahi.com

 東京都渋谷区のNHK放送センター近くの路上で映像制作会社員の男性(48)が首を切りつけられた事件で、警視庁は30日、韓国籍で住所、職業不詳の李宰弦(リジェヒョン)容疑者(46)を殺人未遂の疑いで逮捕し、発表した。李容疑者は事件翌日、「自分がやった」と渋谷署に出頭。その際「無責任な報道をする日本メディアへのメッセージだ」などと話していたというが、逮捕容疑については黙秘しているという。

 捜査1課によると、李容疑者は18日午後9時半ごろ、渋谷区神山町の路上で、同センター内で仕事を終えて出てきた男性の首を刃物で切りつけ、3カ月の重傷を負わせた疑いがある。周辺の防犯カメラの映像で、約4時間前から現場付近をうろつく姿を確認したという。翌19日に出頭した李容疑者は在留資格の期限が切れていたことから、出入国管理法違反(不法残留)容疑で逮捕されていた。昨年12月に福岡県から入国していたという。

渋谷の切りつけ、韓国籍の男逮捕 警視庁、殺人未遂容疑

www.asahi.com

 東京都渋谷区のNHK放送センター近くの路上で映像制作会社員の男性(48)が首を切りつけられた事件で、警視庁は30日、韓国籍で住所、職業不詳の李宰弦(リジェヒョン)容疑者(46)を殺人未遂の疑いで逮捕し、発表した。李容疑者は事件翌日、「自分がやった」と渋谷署に出頭。その際「無責任な報道をする日本メデ…

渋谷の路上切りつけ事件 韓国人の男を殺人未遂容疑で逮捕

5月30日 11時38分

www3.nhk.or.jp

今月、東京・渋谷区の路上で、会社員の男性が刃物で切りつけられて大けがをした事件で、警視庁は、出頭してきた韓国人の46歳の男が事件に関わったとして、殺人未遂の疑いで逮捕しました。

逮捕されたのは、韓国人で、住所と職業がいずれも不詳のリ・ジェヒョン容疑者(李宰弦)(46)です。

今月18日の夜、東京・渋谷区神山町の路上で、神奈川県に住む48歳の会社員の男性が突然、男に後ろから刃物で首を切りつけられて大けがをしました。

事件の翌日にリ容疑者が渋谷警察署に出頭し、事件への関与を認めたため、警視庁が捜査を進めた結果、現場周辺の防犯カメラの映像に写っていたことなどから、30日、殺人未遂の疑いで逮捕しました。

警視庁によりますと、被害者の男性はNHKの業務委託先の映像制作会社の社員で、当時は仕事を終えてNHK放送センターを出た直後だったということです。
リ容疑者は出頭した際の調べに対し、日本のメディアの報道への不満を供述していましたが、逮捕後の調べには「言いたくない。黙秘する」と供述しているということです。

警視庁は事件に至った詳しい経緯や動機についてさらに調べを進めています。

渋谷の路上切りつけ事件 韓国人の男出頭 逮捕へ

5月29日 23時03分

www3.nhk.or.jp

今月、東京・渋谷区の路上で、会社員の男性が刃物で切りつけられて大けがをした事件で、韓国人の46歳の男が警視庁に出頭し、事件への関与を認めていたことが捜査関係者への取材でわかりました。警視庁は男が事件に関わったとみて殺人未遂の疑いで30日にも逮捕する方針です。

今月18日の夜、東京・渋谷区神山町の路上で、神奈川県に住む48歳の会社員の男性が突然、男に後ろから刃物で首を切りつけられて大けがをしました。

警視庁は殺人未遂事件として捜査を進めていたところ、事件の翌日、韓国人の46歳の男が渋谷警察署に出頭し、事件への関与を認めていたことが捜査関係者への取材でわかりました。

警視庁などによりますと、被害者の男性はNHKの業務委託先の映像制作会社の社員で、当時は仕事を終えてNHK放送センターを出た直後だったということです。

警視庁は現場の状況などから男が事件に関わったとみて殺人未遂の疑いで30日にも逮捕する方針です。

渋谷NHK関連社員切りつけ 容疑で韓国籍の男を逮捕 「無責任な報道をするメディアへのメッセージだ」

www.sankei.com

 東京都渋谷区のNHK放送センター近くの路上で今月18日、男性が刃物で切りつけられて首を負傷した事件で、警視庁捜査1課は30日、事件に関与した疑いが強まったとして、殺人未遂容疑で、韓国籍の住所・職業不詳、李宰弦(リ・ジェヒョン)容疑者(46)を逮捕した。調べに対し、黙秘している。

 逮捕容疑は18日午後9時半ごろ、同センターを出た映像制作関連の男性会社員(48)の首に、いきなり刃物で切りつけて殺害しようとしたとしている。男性は首を15センチ切る大けがを負った。

 李容疑者は19日に「自分がやった」と渋谷署に出頭。当初、「無責任な報道をする日本のメディアへのメッセージだ」などと話していたという。警視庁が同日、不法滞在したとして入管難民法違反容疑で逮捕し、捜査していた。

【寝言】バーチャル YouTuber 三千人

vtuber 雨後の筍

戦国四君子の食客三千人。鶏鳴狗盗の士も許されたが。バーチャルYouTuberも三千人!

しかも支那人的な三が沢山、千も沢山という修辞ではなく、マジで三千人越えの模様。近々の一か月で千人参入とかw 参入セットを売った方が儲かる。

www.moguravr.com


Fortran Vtuber も有りか!
fortran66.hatenablog.com


バーチャルYouTuberはじめてみる

バーチャルYouTuberはじめてみる

【朗報】謎の技術で Fortran 内から Python 利用

forpy による Fortran Python 相互運用

Fortran プログラム中から Python のデータ構造を使ったり、モジュールの import やコマンド呼び出しが可能であるとともに、逆に Fortran の shared lib として Python のモジュールを作成できるようです。

コンパイラは gfortran/intel どちらも行けるようです。

github.com

驚異の技術です。詳しくは見ていませんが fypp という fortran preprocessor を利用して Python との Link を可能にする C 言語とのインターフェースを生成しているようです。

github.com

bash on ubuntu on windows 上の intel fortran での実行例

サンプルおよびテストは全て、わずかな修正で動作しました。

Matplotlib の実行

Windows10 上の WSL の bash on ubuntu on windows 上にインストールした intel fortran v.18 によってコンパイル&実行しています。

preprocessor が必要なプログラムになっています。intel fortran の場合、 ifort -fpp xx.f90 と -fpp オプションで乗り切れます。(もしくは拡張子を大文字の F90 にする。)

GitHub - ylikx/forpy: Forpy - use Python from Fortran
f:id:fortran66:20180522014342p:plain

intrinsic modules

intel fortranコンパイルするときのオプションは、例題に示されている gfortran のものと同じで良いようです。ただし、ソースプログラム中の、use iso_fortran_env や use iso_c_binding を、use, intrinsic :: iso_fortran_env や use, intrinsic :: iso_c_binding に直す必要がありました。

なお、ファイル拡張子が F90 と大文字になっているので f90 の小文字派の人はうっかりに注意ですw
[追記] 大文字の F90 拡張子は、プリプロセッサを稼働させるそうです。

tests

また tests フォルダー中の Makefike では、コンパイラの指定を ifort にしました。上述のエラーの出るソースファイルでは intrinsic module の修正をする必要があります。 make runtests で全テストを合格しました。

python module

https://github.com/ylikx/forpy#developing-python-modules-in-fortran

元のサンプルのままでは segmentation fault を起こしますが、関数 print_args の宣言行の末尾に bind(c, name = 'print_args') をつけてやれば動くようになります。

function print_args(self_ptr, args_ptr, kwargs_ptr) result(r) bind(c, name = 'print_args')
O@HP8:~/forpy/DLL$ ifort -c -fpic ../forpy_mod.F90
O@HP8:~/forpy/DLL$ ifort -shared -fpic -o extexample01.so extexample01.f90 forpy_mod.o `python3-config --ldflags`
O@HP8:~/forpy/DLL$ python extexample01.py
Segmentation fault (core dumped)
O@HP8:~/forpy/DLL$ vi extexample01.f90
O@HP8:~/forpy/DLL$ ifort -shared -fpic -o extexample01.so extexample01.f90 forpy_mod.o `python3-config --ldflags`
O@HP8:~/forpy/DLL$ python extexample01.py
('hello', 42)
{'key': 'abc'}
3.141592653589793
O@HP8:~/forpy/DLL$

gfortran ・・・うまくゆかず

なお gfortran Ver.7 ではコンパイル時にエラーが出て実行できませんでした。

O@HP8:~/forpy/gfort5$ ls
forpy_mod.F90  forpy_mod.mod  forpy_mod.o
O@HP8:~/forpy/gfort5$ gfortran -c -fPIC ../forpy_mod.F90
O@HP8:~/forpy/gfort5$ gfortran ../intro_to_forpy.F90 forpy_mod.o `python3-config --ldflags`
lto1: internal compiler error: in lto_tag_to_tree_code, at lto-streamer.h:1005
Please submit a full bug report,
with preprocessed source if appropriate.
See <file:///usr/share/doc/gcc-7/README.Bugs> for instructions.
lto-wrapper: fatal error: gfortran returned 1 exit status
compilation terminated.
/usr/bin/x86_64-linux-gnu-ld: error: lto-wrapper failed
collect2: error: ld returned 1 exit status
O@HP8:~/forpy/gfort5$ gfortran --version
GNU Fortran (Ubuntu 7.3.0-16ubuntu3~16.04.1) 7.3.0
Copyright (C) 2017 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

O@HP8:~/forpy/gfort5$

追記 (H30.8.17)

gfortran でうまくいかなかったのは、Python を anaconda で入れた場合の症状のようです。対処法が出ていました。たしかにうまくゆきます。

GitHub - ylikx/forpy: Forpy - use Python from Fortran

Using forpy with Anaconda
When using forpy with Anaconda and gfortran, you might encounter the following error:

/usr/bin/x86_64-linux-gnu-ld: error: lto-wrapper failed
collect2: error: ld returned 1 exit status
A solution to this problem is to add the -fno-lto compiler flag in the linking step:

gfortran -c forpy_mod.F90
gfortran intro_to_forpy.F90 forpy_mod.o -fno-lto `python3-config --ldflags`

蛇はマジむかつくわw

邪神ちゃんドロップキック 邪神ちゃん 全高約135mm PVC製 塗装済み 完成品 フィギュア

邪神ちゃんドロップキック 邪神ちゃん 全高約135mm PVC製 塗装済み 完成品 フィギュア