プロフィールPROFILE

裏 RjpWikiさんのプロフィール

住所
未設定
出身
未設定

自由文未設定

ブログタイトル
裏 RjpWiki
ブログURL
https://blog.goo.ne.jp/r-de-r
ブログ紹介文
文字通り,RjpWiki の裏を行きます
更新頻度(1年)

55回 / 365日(平均1.1回/週)

ブログ村参加:2014/11/25

裏 RjpWikiさんのプロフィール
読者になる

裏 RjpWikiさんの人気ランキング

  • IN
  • OUT
  • PV
今日 11/15 11/14 11/13 11/12 11/11 11/10 全参加数
総合ランキング(IN) 圏外 圏外 圏外 圏外 圏外 圏外 圏外 978,010サイト
INポイント 0 0 0 0 0 0 0 0/週
OUTポイント 0 0 0 0 0 0 0 0/週
PVポイント 0 0 0 0 0 0 0 0/週
科学ブログ 圏外 圏外 圏外 圏外 圏外 圏外 圏外 2,265サイト
数学 圏外 圏外 圏外 圏外 圏外 圏外 圏外 161サイト
今日 11/15 11/14 11/13 11/12 11/11 11/10 全参加数
総合ランキング(OUT) 圏外 圏外 圏外 圏外 圏外 圏外 圏外 978,010サイト
INポイント 0 0 0 0 0 0 0 0/週
OUTポイント 0 0 0 0 0 0 0 0/週
PVポイント 0 0 0 0 0 0 0 0/週
科学ブログ 圏外 圏外 圏外 圏外 圏外 圏外 圏外 2,265サイト
数学 圏外 圏外 圏外 圏外 圏外 圏外 圏外 161サイト
今日 11/15 11/14 11/13 11/12 11/11 11/10 全参加数
総合ランキング(PV) 圏外 圏外 圏外 圏外 圏外 圏外 圏外 978,010サイト
INポイント 0 0 0 0 0 0 0 0/週
OUTポイント 0 0 0 0 0 0 0 0/週
PVポイント 0 0 0 0 0 0 0 0/週
科学ブログ 圏外 圏外 圏外 圏外 圏外 圏外 圏外 2,265サイト
数学 圏外 圏外 圏外 圏外 圏外 圏外 圏外 161サイト

新機能の「ブログリーダー」を活用して、裏 RjpWikiさんの読者になりませんか?

ハンドル名
裏 RjpWikiさん
ブログタイトル
裏 RjpWiki
更新頻度
55回 / 365日(平均1.1回/週)
読者になる
裏 RjpWiki

裏 RjpWikiさんの新着記事

1件〜30件

  • 消費増税ーーなんで,理性的な判断ができないか?

    今回の消費増税にあたって,駆け込み購入が目立たないという報道が多数あったが,ここ数日,それを煽るような報道が多く見られるようになった。たとえば,リクルートスーツ。定価2万5千円を駆け込み購入するとか。従来の消費税なら税率8%で,税込み2万7千円。消費税分は2000円。増税後は税率10%で,税込み2万7千500円。消費税分は2500円。差額は500円。これは大きいか?リクルートスーツ。就職活動の後も使うだろう。何年間?まあ,2年としても,2年間の内,500円くらいのラーメンを一回我慢するだけで500円はカバーできる。2年間の内の一回のラーメン代も我慢できないというなら,それは仕方ない。普通の人ならできるだろう。駆け込み購入は,トイレットペーパーから猫の餌までということだ。前回の消費増税のとき,ふりかけ(ごはんにか...消費増税ーーなんで,理性的な判断ができないか?

  • 多倍長精度による円周率 π の計算

    「多倍長精度による計算」に,Rによる円周率πの計算プログラムが掲示されている。Rでもgmpパッケージを使えば多倍長制度による計算はできるので,このプログラムを書き換えればよい。Pythonでもgmpyモジュールを使えば多倍長精度による計算ができるが,例によって64ビット整数型を使えば,簡単にプログラムを移植できる。ということで,以下のようなプログラムになる。移植に際して注意すべきことは,除算は/ではなく,整数除算//を使うこと。また,πは浮動小数点すうであるが,指定した桁の固定小数点数として扱う。具体的には小数点以下の桁数(acculacyで指定)の下に小数点を置いて,整数値πxacculacyを計算する。defcalc_pi(acculacy=10000):FACTOR=10**acculacypi=0n,a...多倍長精度による円周率πの計算

  • 長精度整数演算

    Sumofthreecubesfor42finallysolved--usingreallifeplanetarycomputer3つの整数の3乗和が42になる解が見つかったとのこと。なんで,42なのかは興味がないので説明しない。かなり大きな整数なので,Rの普通の変数型numeric(double)では扱えないので,gmpライブラリを使う必要がある。>library(gmp)>x=as.bigz("-80538738812075974")>library(gmp)>x=as.bigz("-80538738812075974")>y=as.bigz("80435758145817515")>z=as.bigz("12602123297335631")>ans=x^3+y^3+z^3>ansBigInteger('...長精度整数演算

  • マシュマロを食べることができるだろうか?

    本日(2019/09/08)の日本テレビ,世界の果てまで行ってQで,100mの高さから10個くらいのマシュマロを一度に投げて,下で待ち受ける8人ほどが口に入れることができるだろうかという実験をやっていた。結果は失敗だったけど,どんな条件だと成功する確率がどれくらいになるのかシミュレーションしてみた。library(MASS)#9人を40センチ間隔の格子点に配置(単位はセンチメートル)x=c(-40,-40,-40,0,0,0,40,40,40)y=c(-40,0,40,-40,0,40,-40,0,40)person=length(x)sigma=4000#二次元正規分布の分散marshmallow=10#一度に投げるマシュマロの数trial=1000#実験回数result=integer(trial)#マシュ...マシュマロを食べることができるだろうか?

  • Python って奴は,なんて自由な奴なんだ!

    プログラミングの初歩の段階で,「1より大きくて10より小さい場合は'OK',そうでない場合は'NO'を書くぷろぐらむを作りなさい」などというのがあり,初心者の中には,if(1<x<10)print('OK')などと書きたくなる人がいる。大抵の言語はこれはだめで,if(1<x&&x<10)print('OK')としないといけない。しかし,である。Pythonはこれを許す。なんですと!!ちっとも知らなかった。defif_test(x):if1<x<10:print('OK')else:print('NG')if_test(0)#NGif_test(5)#OKif_test(10)#NGPythonって奴は,なんて自由な奴なんだ!

  • ビットカウント(ビットが 1 である個数)

    「ビットカウントする高速アルゴリズムをPythonで実装しながら詳しく解説してみる」というのを見つけた。Pythonの整数は64ビット整数,Rは32ビット整数であることを考慮し,以下のように翻訳した。bitcount=function(n){#nの立っているビット数をカウントする関数(nは64bit整数)#2bitごとの組に分け、立っているビット数を2bitで表現するn=n-bitwAnd(bitwShiftR(n,1),0x55555555)#4bit整数に上位2bit+下位2bitを計算した値を入れるn=bitwAnd(n,0x33333333)+bitwAnd(bitwShiftR(n,2),0x33333333)n=bitwAnd((n+bitwShiftR(n,4)),0x0f0f0f0f)#8bit...ビットカウント(ビットが1である個数)

  • utils:::combn の代替関数 next.combn

    utils:::combnは組み合わせを配列で返すが,場合によっては一つずつ返してくれる方がうれしいこともある。そこで,以下に示すnext.combnを書いてみた。next_combn=function(n,r,a){t=rwhile(t>=1&&a[t]==n-r+t){t=t-1}if(t==0){return(FALSE)}a[t]=a[t]+1for(iint:r){a[i]=a[t]+i-t}return(a)}使い方は,n,rを指定してaを1,2,...,nのベクトルに定義してn=5r=3a=1:n以下のように引用する。print(a[1:r])の所に,本来意図している処理プログラムを書く。while(TRUE){print(a[1:r])a=next_combn(n,r,a)if(length(a...utils:::combnの代替関数next.combn

  • e1071:::permutations の代替関数 next.permutation

    e1071ライブラリにはpermutations関数がある。permutationsは処理が速いが,配列で返すために大きなメモリーを消費する。時間はかかってもよいので,メモリー消費のない関数が必要なこともある。そのようなときのために以下の関数を書いた。next.permutation=function(x){n=length(x)pos1=0for(iin(n-1):1){if(x[i]<x[i+1]){pos1=ibreak}}if(pos1==0){return(0)}for(iinn:(pos1+1)){if(x[pos1]<x[i]){pos2=ibreak}}temp=x[pos1]x[pos1]=x[pos2]x[pos2]=templen=(n-pos1)%/%2if(len>=1){for(ii...e1071:::permutationsの代替関数next.permutation

  • bincombinations から n-combinations まで

    Rのe1071ライブラリにbincombinationsがある。bincombinations(3)で,2^3×3行列の結果が返される。>bincombinations(3)[,1][,2][,3][1,]000[2,]001[3,]010[4,]011[5,]100[6,]101[7,]110[8,]111要するに,n桁の2進数の各桁が2^n行で表される。ただ,これだと結果が行列で返されるので,不都合な場合もある。また,一般にn進数の結果が必要になることもある。bincombinationsを一般化するのは,このページの下の方に再掲載しておく。まず,与えられたn進数表記の次の数の表記を与える関数を書いてみる。next.Ncombination=function(x,base=2){n=length(x)car...bincombinationsからn-combinationsまで

  • insert, pop, count, index, concat, remove を定義する

    Rには,ベクトルに要素を挿入するappendがあるが,関連する他の関数がないので作ってみた。1.要素の挿入index(x,a,pos=length(x))ベクトルxのpos位置にaを挿入する。先頭に加える場合はpos=0とする。insert=function(x,a,pos=length(x)){assign(deparse(substitute(x)),append(x,a,after=pos),envir=.GlobalEnv)}>x=c('w','q','p','w')>insert(x,'first',0)>insert(x,'last')>x[1]"first""w""q""p""w""last"2.要素の取り出しpop(x,pos=length(x))ベクトルxのpos位置にある要素を取除く。取り除...insert,pop,count,index,concat,removeを定義する

  • 連続変数をカテゴリー化する

    連続変数をカテゴリー化する関数としてはfindInterval関数がある。set.seed(123)x=rnorm(100000)v=c(-Inf,-1.2,-0.8,0,0.5,1.4,2.2,Inf)system.time({a=findInterval(x,v)table(a)})もし,以下のような関数を作ってみると,findInterval関数より,15倍遅い。system.time({b=sapply(x,function(y)sum(y>v))table(b)})連続変数をカテゴリー化する

  • 四捨五入関数は round ではない(かも)

    三省堂大辞林によれば,四捨五入とは,「必要とする位の次の位の数字が四以下ならばこれを切り捨て、五以上ならば切り上げて、上の位に一を加える方法」とある。Rでは,round関数がこれを行う。1.5を,小数第1位で四捨五入を行うと2である。>round(1.5)[1]2しかし0.5を,小数第1位で四捨五入を行うと0である。>round(0.5)[1]0よって,roundは四捨五入を行う関数ではない。roundは,偶数への丸め(最近接丸め,JIS丸め,ISO丸め,銀行丸めなどとも呼ばれる)を行う関数である。これは,「端数が0.5より小さいなら切り捨て、端数が0.5より大きいならは切り上げ、端数がちょうど0.5なら切り捨てと切り上げのうち結果が偶数となる方へ丸める。(JISZ8401の規則A)」として定められており、規則...四捨五入関数はroundではない(かも)

  • なかなか頷きがたい

    数学的に考える「じゃんけん」で有利な手は何か>725人から集めたのべ11567回のじゃんけんデータからみたじゃんけんの確率はどうだったのか?集計結果は以下である。グーが4054回、チョキが3664回、パーが3849回。したがって、「グーを出す人が多いことから一般にじゃんけんではパーが有利」だと導かれた。それぞれ,35.0%,31.7%,33.3%である。一番弱いチョキと比べてグーだと3%しか有利ではない。確かに,のべ11567回の結果からはp値=5.166e-05などと言うことになるんだが,これはいかに何でもサンプルサイズが大きすぎる。それぞれの手の割合が同じでも,サンプルサイズを1/4にしたら,p値=0.0826で有意ではない。サンプルサイズが2892でも相当大きい。ねえ!世論調査や視聴率調査のときのサンプル...なかなか頷きがたい

  • Mac でのファイル名(UTF-8-MAC)

    こっちにも控えておこうlist.files()で読み込んだファイル名に濁音,半濁音が含まれるときの対処法>y<-list.files()>x=c('織田信長','パンダ')>x%in%y[1]TRUEFALSE>x=iconv(x,"UTF-8","UTF-8-MAC")>x%in%y[1]TRUETRUEMacでのファイル名(UTF-8-MAC)

  • 二次方程式の解を求めるプログラムのレビュー

    二次方程式の解を求めるプログラム例はネット上でもよく見られる。本来(?)お手本を示すべきものであるが,問題を含むものが少なくない(多いということだ)。結論を先に述べておく。(0)全部で40サイトを対象とした(1)floatを使っているのは8サイト。そのうちCによるものは6サイト。Cは12サイトだが,そのうちfloatを使っているのは6サイト。Cはよく使われているが,不適切なものも多い。(2)桁落ちを考慮しているサイトは8サイト(20%)。(3)虚数型を使用しているのは3サイトに過ぎない。(4)虚数解の場合を対象としていないのは17サイト。実部虚部を表示するのは19サイト。==以下は各サイトのレビュー==●●https://webkaru.net/clang/quadratic-formula/はC言語によるもの...二次方程式の解を求めるプログラムのレビュー

  • 紙は何回折れるか?

    かなり有名な事実で,「ほぼ8回」というのが正解。厚さ0.1mm,一辺100mの紙を,偶数回ごとに正方形になるように8回折れば,厚さは0.1/1000*2^8=0.0256mで,正方形の一辺の長さは6.25m。まだ十分に折れそうではあるが,10回折ると,厚さは0.1024m,一辺の長さ3.125mだが,折られた紙のそれまでの折り線付近が折りにくくなる。フットボール競技場大の紙を折るときは11回折れたという動画がある。柔らかい紙を一方向に折っていくともう少し多く折れる。1200メートルのトイレットペーパーを半分ずつに折るのは12回できたという写真がある(ワンカットの動画でないので正しさの証明はないが)。このとき,最終のトイレットペーパーの長さは1200*0.5^12=0.2929688で,まあ写真にあるように肘の長...紙は何回折れるか?

  • 2直線の交点座標を求めるプログラム(Python 版)

    何故か相も変わらず「2直線の交点座標を求めるプログラム」がよく参照されているそこにあるのはRプログラムなので,Pythonで書いたものをおまけに載せておこうdefprog(x1,y1,x2,y2,x3,y3,x4,y4):a1=(y2-y1)/(x2-x1)a3=(y4-y3)/(x4-x3)x=(a1*x1-y1-a3*x3+y3)/(a1-a3)returnx,(y2-y1)/(x2-x1)*(x-x1)+y1prog(4,8,8,0,2,2,12,7)なんでこんなものが検索されるのやら?2直線の交点座標を求めるプログラム(Python版)

  • Cohen の d

    中澤さんが>effsizeパッケージのcohen.d()関数の計算は何か変だ。大久保・岡田(2012)のSpの式で標本分散のはずのところが,不偏分散になっている。どちらかがCohenの定義を誤解していると思われるが,不偏分散にサンプルサイズを掛けるのは筋が通らないので,たぶんeffsizeパッケージが間違っているのだろう。と書いている。大久保・岡田(2012)がないのでよく分からないのだが,Cohen'sdは,奥村先生のページにもあるが,で,stdev=sqrt(((n1-1)*s[1]^2+(n2-1)*s[2]^2)/(n1+n2-2))だからいいのではないか?それとも,私が何か勘違いしている?>x=c(1,2,3,2,3,4)>y=c(5,4,3,4,5,6,5)>cohen.d(c(x,y),facto...Cohenのd

  • ぼったくり宝くじ

    オイ,録音してないやろななんですか,いきなり。してませんがな。なんでせなならんの。ちょっとヤバい話だ。良く効け。おれらも相当あくどい仕事してるけど,世の中には堂々とひどい商売やってる奴がいるぜ。なんですか?宝くじや。あ,今,「買わないという選択はない」とかコマーシャルで言ってる。そや。でも兄貴。あのお堅い銀行さんも噛んでいる,まともな商売じゃないんですか。どこがまともじゃ。アホンダラ。ええか,100万円分の宝くじ買ったとするやろ,なんぼ当たるかわかるか?え〜っ。100万円も買うんでっか。1000万位は当たりそう。あほ抜かせ。ええとこ行って,46万位しか当たらんのじゃ。え〜っ。え〜っ。なんでそんな。ええか。ちゃんと,総務省,総務省やぞ,のWebページというかPDFhttp://www.soumu.go.jp/ma...ぼったくり宝くじ

  • 貯蓄額倍増?

    ご隠居さん,コンチハ。ああ,はっつあんか。まあお上がり。ちょっとご隠居さんにうかがいたいことがあるんですけど。なんでも聞いとくれ。「定年後に2000万円必要だ!」なんて言ってる人がいるんだけど,ホントですか。まあ,ほんとのようなうそのような。なんだろね。今,銀行の利率なんて雀も気の毒がるほど低いじゃないですか。そうだねえ,昔は年利率5%なんてときもあったんだが,今じゃ0.2%なんてことも普通だねえ。あっしが,今100万円持ってるとして,年利0.2%の複利で預けるとして,元利合計が200万円になるのは何年後です?おや,驚くねえ。100万も持ってるのかい?冗談言っちゃいけね〜。持ってるとしてといったでしょうに。そうさなあ,年利0.2%なら350年ぐらい掛かるな。ご冗談でしょ。その前にテメーの命がなくならあ。そうさな...貯蓄額倍増?

  • 土用ウナギの日

    土用の丑の日に限らないが,よく食の老舗で「江戸時代から継ぎ足し継ぎ足ししてきたタレ」とかいう話があり,そのたびに,「ばかか」と思うまあ,ウナギのタレがどういう単位で存在するのかはさておき,たとえば単位不明でn=1000000としよう。毎日毎日そのうちの千分の1であるm=1000を使い,そのぶんだけ新しいものを補充するとしよう。当初,n個の要素は全て1を持っている。a=rep(1,n)毎日毎日そのうちのr=sample(n,m,replace=TRUE)の要素が新しいものに置き換えられる。つまり,使われた要素は1から0に変えられる。a[r]=0一日が終わったとき,当初の要素が幾つあるかはsum(a)でわかる。というようなシミュレーションプログラムが以下のようになるであろう。n=1000000m=1000a=rep...土用ウナギの日

  • トリボナッチ数列

    「再帰関数と非再帰関数」,「事象が連続発生するまでの試行回数」にも出てきたトリボナッチ数列だが,「n階段があるとき,1段昇る,2段昇る,3段昇るのどれかを選択して昇るとき何通りの昇り方があるか」という設問がされることもある。よく,「アルゴリズムをわかりやすく表現できる再帰プログラム」の例としても出てくるのが以下のプログラム。フィボナッチ数列の拡張ということだ。defstep(n):ifn==0:return1ifn<0:return0returnstep(n-1)+step(n-2)+step(n-3)でも,これは相当スタックも喰うし,実行時間も馬鹿にならない。importtimes=time.time()foriinrange(30):print(i,step(i))print(time.time()-s)3...トリボナッチ数列

  • 事象が連続発生するまでの試行回数

    やあ,ジュン。元気かい?元気だよ。前に,ポケモンGoのことで教えてもらったけど,またいいかな?いいよ。でも,ポケモンGoに興味がない人にも分かるような例で説明してくれる?そうか〜。うんとね。あ,そうそう。野球にたとえるとこんなになるかな。「3割バッターが,3連続ヒットを打てるまでに何回打席に立たなければならないか」でいいかな。前に教えてもらった負の二項分布に似ている気もする。ヒットを打ったら1,打たなかったら0で表すと,打席の記録は100110100111のように最後の3桁が111になったら3連続ヒットを打ったということ。当然だけど,最後の3桁より左には1が3つ以上連続することはない。この例だと,12打席目で3連続ヒットを達成ということだけど,確率変数xは3連続ヒットまでの打席数として9と考えるほうがいいんだろ...事象が連続発生するまでの試行回数

  • デカルトの正葉線

    デカルトの正葉線foliumofDescartesを描いてみようと思い,やってみたが自力ではちょっと難しかった。媒介変数表示だと,大抵はx=3*a*t/(1+t^3)y=3*a*t^2/(1+t^3)が紹介されている。しかしこれだと,原点付近が描画できない。そこで,更にt=(1+s)/(1-s)と変数変換して,x=1.5*a*(1-s)*(1-s^2)/(1+3*s^2)y=1.5*a*(1+s)*(1-s^2)/(1+3*s^2)で描画する。(ということだ)s=seq(-3.5,3.5,by=0.01)x=1.5*a*(1-s)*(1-s^2)/(1+3*s^2)y=1.5*a*(1+s)*(1-s^2)/(1+3*s^2)plot(x,y,type="l",asp=1)abline(-a,-1,col=2,...デカルトの正葉線

  • ベアストウ法による求解

    別名ヒチッコック法(というか,ヒチコック・ベアストウ法として覚えていたが)により,実係数の高次代数方程式を解く。虚数解になる場合も含み,n次式の場合はn個の解を全て求める。インターネット上にもプログラム例が見られるが,数値計算的,プログラミング的に見た場合,ちょっとひどいなあというのも(箇所も)見られるので,お手本とはいわないまでも,まずまずのプログラムを示しておこう。まずはRによるプログラム例,その後にPythonによるプログラム例を掲示する。quadratic.root=function(p,q){d=p^2-4*qif(d>=0){x1=(-p-ifelse(p>=0,1,-1)*sqrt(d))/2#絶対値の大きいほうの解を求めるx2=q/x1#もう一方は解と係数の関係から得る}else{x1=(-p+...ベアストウ法による求解

  • グラフの軸の目盛り設定

    中澤さんが>出勤前にふと思いついて,trux.Rというコードを書いてみた。これで問題なければ,人口ピラミッド描画パッケージのデフォルト横軸目盛設定アルゴリズムはこれに変えようと思う。trux<-function(X){SC<-10^as.integer(log10(X)-1)MX<-(X%/%(6*SC))*SC*2return(0:4*MX)}trux(132)trux(495)trux(8734)trux(13044050)trux(94822324)>trux(132)[1]04080120160>trux(495)[1]0160320480640>trux(8734)[1]028005600840011200>trux(13044050)[1]0.0e+004.0e+068.0e+061.2e+071....グラフの軸の目盛り設定

  • ポケモン Go と負の二項分布

    よう!ヨッシーじゃないか。歩きながらポケモンやってると危ないよ。おう!ジュンか。ついつい夢中になっちゃったんだけど,反省しないといけないね。なんか,おもしろい話はないかい?そうだねえ。ポケモンといえば,最近ポケモンの写真を何回か撮ると,何回目かに特別なポケモンが写真に写り込んでその特別なポケモンをゲットできるという機能が実装されたんだけど,1回でゲットできることはほとんどなくて,大抵2,30回目にやっとゲットできるようなんだ。ゲットできるまでに何回くらい写真を撮らないといけないかは,1回あたりのポケモンゲットの確率とどんな関係にあるんだろうね?「確率pで起きる事象がk回起きるまでに,その事象が何回起きなかったか」という分布を負の二項分布というというのがあったね。負の二項分布か〜。その分布は何で決まるの?pとkと...ポケモンGoと負の二項分布

  • アルキメデスの方法による円周率の近似値

    直径1の円の外接正n角形の周長をPn,内接正n角形の周長をpnとすれば,n=3*2^k,k=0,1,2,...P2n=2*Pn*pnp2n=sqrt(P2n*pn)ということで,以下のプログラムimportscipyasspdefPi_Archimedes():print("{0:>2s}{1:>9s}{2:>17s}{3:>17s}".format("k","n","Pn","pn"))k=0#正三角形から始めるPn=3*sp.sqrt(3)pn=Pn/2whileTrue:n=3*2**kprint("{0:2d}{1:9d}{2:.15f}{3:.15f}".format(k,n,Pn,pn))ifPn==pn:returnPnk+=1P2n=2*Pn*pn/(Pn+pn)p2n=sp.sqrt(P2n*...アルキメデスの方法による円周率の近似値

  • ロンバーグ積分(その2)

    前報の「ロンバーグ積分」には,Cで書かれたプログラムをPythonに書き換えたものを示した。両言語は共に「0オリジン」https://jp.mathworks.com/matlabcentral/fileexchange/34-rombint-mには,MATLABで書かれたプログラムがある。#ROMBINTNumericalevaluationofanintegralusingtheRombergmethod.##Q=rombint(F,A,B)approximatestheintegralofF(X)fromAtoBto#withinarelativeerrorof1e-10usingRomberg'smethodofintegration.#Fisafunction.Thefunctionmustretur...ロンバーグ積分(その2)

  • ロンバーグ積分

    ロンバーグ積分とは,まず合成台形公式で近似し、次にRichardsonの補外法で近似をより正確にする。Romberg'smethodに,Cで書かれたプログラムがある。これを,Pythonに書き換えると以下のようになる。forをベクトル計算にするなどにより,ずいぶん短くなる?importscipyasspdefRomberg(f,a,b,max_steps=1000,acc=1e-14):Rp=sp.empty(max_steps)#previousrowRc=sp.empty(max_steps)#currentrowh=b-aRp[0]=(f(a)+f(b))*h/2foriinrange(1,max_steps):h/=2Rc[0]=h*sum(f(sp.arange(1,2**i,2)*h))+Rp[0]...ロンバーグ積分

カテゴリー一覧
商用