算額あれこれ

SymPy で二重根号を外す

以前の記事で，「SymPyでは二重根号を外せない」と書いたが，そんなことはなかった。sqrtdenest()があった。Maximaにもあるようだが，SymPyのsqrtdenest()に言及した日本語の記事は1つしか見つからなかった。ちなみに，そこには「二重根号外しを二重にやらなければいけないケースにどのくらい対応できているか試してみた」とあり，うまく行かないという以下の例が示されていた。Pythonの場合importsympyasSS.init_printing()a=S.sqrt(S.sqrt(49+20*S.sqrt(6)))aS.sqrtdenest(a)うまくやるためにはsqrtdenestを2回使うことである（Pythonのsqrtdenestのソースドキュメントにヒントが書いてあった）。S....SymPyで二重根号を外す

2022/09/30 17:41

大学入試問題を SymPy で解く

以下にある問題をSymPyで解いてみる。大学入試数学の問題http://mikiotaniguchi.com/main/center/main.htm22092501剰余の定理22京都産業大22082401定積分01芝浦工業大14121806円と直線05大阪市立大大学入試問題をSymPyで解く

2022/09/27 15:19

バーゼル問題

1/n^2をn=1,2,...,∞として，和を取るとπ^2/6に収束する。usingSymPysummation(1/n^2,(n,1,oo)) >string"pi^2/6"バーゼル問題

2022/09/26 15:10

ゲルフォントの定数

julia>ℯ^π23.140692632779267julia>big(ℯ)^π23.14069263277926900572908636794854738026610624260021199344504640952434235069032有効桁をほぼ300桁になるように設定する。julia>setprecision(BigFloat,1000)1000julia>big(ℯ)^π23.14069263277926900572908636794854738026610624260021199344504640952434235069045278351697199706754921967595270480108777314442804441469383584471744587960984936532796...ゲルフォントの定数

2022/09/26 13:20

年齢の計算

年齢の計算生年月日year,month,dayとして，ある時点year2,month2,day2での年齢を求めるのに以下のようなプログラムを書いていた。関数定義age(year,month,day,year2,month2,day2)=year2-year-(month>month2 (month==month2&&day>day2))実行age(2020,9,19,2022,8,29),age(2020,9,19,2022,9,18),age(2020,9,19,2022,10,18),age(2020,9,19,2022,9,19)実行結果(1,1,2,2)日付をyyyymmddのような8桁の整数で受け渡しすると，以下のように簡単なプログラムになる。ちゃんと上のプログラムと同じように月と日の大小順を...年齢の計算

2022/09/19 18:06

Julia で閏年判定

Datesパッケージのisleapyear()を使うjulia>usingDatesjulia>isleapyear.([2020,2022,2100,2400])4-elementBitVector:1001Juliaで閏年判定

2022/09/19 11:40

Rump の例題--その2 Julia の SymPy で解く

JuliaのSymPyは，Pythonのsympyを利用している。Pythonの整数型は必要に応じて長精度整数になる。求める式には実数が含まれるが，この式を4倍することで整式になる。Rumpの例題--その2JuliaのSymPyで解く

2022/09/13 11:05

Rump の例題

Rumpの例題https://www.tuhh.de/ti3/paper/rump/Ru88a.pdf以下の計算を行うRf=function(a,b)333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/(2*b)library(gmp)a=as.bigq(77617)b=as.bigq(33096)ans=f(a,b)as.numeric(ans)#-0.8273960599468214Juliaf(a,b)=333.75b^6+a^2*(11a^2*b^2-b^6-121b^4-2)+5.5b^8+a/2ba=big(77617)b=big(33096)f(a,b)#-0.8273960599468213681411650954798162919990...Rumpの例題

2022/09/12 16:05

Julia 1.8.1 が公開されました

Julia1.8.1が公開されましたVersion1.8.1(2022-09-06)Officialhttps://julialang.org/releaseJulia1.8.1が公開されました

2022/09/11 21:33

AWK と R のコラボレーション 定式化された入力ファイルからデータファイルへ

一般解はhttps://stackoverflow.com/questions/69385888/parse-multi-line-result-text-file-with-awk-bash-or-rにあった。AWKスクリプトを平易に書き変えた。★以下のような定式化された入力ファイルmelting.txtTheMELTINGresultsare:Enthalpy:-181,100cal/mol(-756,998J/mol)Entropy:-467.3cal/mol-K(-1,953.31J/mol-K)Meltingtemperature:75.13degreesC.TheMELTINGresultsare:Enthalpy:-170,800cal/mol(-713,944J/mol)Entropy:-...AWKとRのコラボレーション定式化された入力ファイルからデータファイルへ

2022/09/07 17:13

AWK で製表

以下のような入力ファイルCodeFormatCodeFormatExample0dd-mmm-yyyyHH:MM:SS07-Sep-200015:38:091dd-mmm-yyyy07-Sep-20002mm/dd/yy09/07/003mmmSep4mSからデータを読み込みCode Format Example:-- :-- :--0 dd-mmm-yyyyHH:MM:SS 07-Sep-200015:38:091 dd-mmm-yyyy 07-Sep-20002 mm/dd/yy 09/07/003 mmm Sep4 m Sのように整形して出力ファイルtableに書き込む。つまり，最終形の表は3列からなり，各列は で区切られている。元のデータは3行で表の1行になる。awkで書くとこんなスクリプト{ge...AWKで製表

2022/09/07 14:42

RStudio で awk を使う

前の記事で```{awk}...```が動かないと書いた。https://gedevan-aleksizde.github.io/rmarkdown-cookbook/custom-engine.htmlの「15.1カスタム言語エンジンを登録する」を読んで，以下のように解決した（のかな?）。最初に一回だけ，以下のRスクリプトを実行する。knitr::knit_engines$set(awk=function(options){code<-paste(options$code,collapse='\n')out<-system2('awk',c('-c',shQuote(code)),stdout=TRUE)knitr::engine_output(options,code,out)})この後，awkによる...RStudioでawkを使う

2022/09/03 15:31

Rstudio で Julia を使う

RStudioでJuliaを使う方法を述べる。とは言っても，RStudioのノートブック作成機能を使うだけである。通常はR言語を選択するために{r}となっているところを{julia}にするだけだ。1.まず，RStudioのファイルメニューから，File-->NewFile-->RNotebookでノートブックウインドウを出す。2.以下の例のように，コードブロック```{julia}...```で挟んでJuliaプログラムを書く。3.右上にあるRunボタンをクリックし，RunCurrentChunkをクリックする。実行すべきチャンクを選択してshift+command+returnでも良い。4.実行結果はチャンクの下に書かれる。{}内に書けるものは，R,Python,Julia...RstudioでJuliaを使う

2022/09/02 09:30

Octave でトーラスを描いてみた

Octaveでトーラスを描いてみた。Juliaに翻訳--219トーラス，3次元プロット参照Juliaのrange(),repeat()に相応する関数がないので，自前の関数を書く。function[rng]=range_(begin_,end_,n)step=(end_-begin_)/(n-1);rng=begin_:step:end_;endfunction[res]=repeat(x,n)res=x;fori=2:nres=[res,x];endendfunctiontorus(n=100,a=2,b=9)u=reshape(repeat(range_(0,2*pi,length=n),n),n,n);v=u';x=(b+a*cos(u))*cos(v);y=(b+a*cos(u))*sin(v);z...Octaveでトーラスを描いてみた

2022/09/01 16:44