einsum!を作る
この記事は数値計算Advent Calendar 2022の1日目の記事です。 前回の記事では既存実装としてNumPyにおけるnumpy.einsumの仕様を見ていきましたが、今回はRustのndarray crate向けにeinsumを作っていきます。完成形は次のようになります： use ndarray::array; use einsum_derive::einsum; let a = array![ [1.0, 2.0], [3.0, 4.0] ]; let b = array![ [1.0, 2.0], [3.0, 4.0] ]; let c = einsum!("ij,jk->ik", a, b); assert_eq!(c, array![ [6.0, 8.0], [12.0, 16.0] ]); このeinsum_derive crateは下記のリポジトリで開発されています：
p進Newton法で平方数判定 - Qiita
浮動小数点数演算を使わずに高速に平方数を判定します。 (実装例はRustですがC++でも書けると思います) use std::num::Wrapping; pub fn is_square(x: u64) -> bool { let xtz = x.trailing_zeros(); if x == 0 { return true; } else if xtz & 1 != 0 { return false; } let x = Wrapping(x >> xtz); let mut y = Wrapping(1_u64); for _ in 0..5 { y = ((Wrapping(3_u64) - y * y * x) * y) >> 1; } let mut xrt = x * y; if xrt.0 & (1 << 32) != 0 { xrt = -xrt; } let xrt
Implementing cosine in C from scratch
Implementing cosine in C from scratch 7/19/2020 Update 7/20: See the discussion of this post on Reddit. Update 3/22: See more discussion of this post on Hacker News and Reddit. Update 6/23: See even more discussion of this post on Hacker News. TL;DR: I explored how to implement cosine using a few different approaches. One of the implementations is nearly 3x faster than math.h if you are ok with 4
非圧縮性 Navier-Stokes 方程式の数値解法２：拡散項の陰解法 - Qiita
はじめに Navier-Stokes 方程式を数値的に解くシリーズ、第２回目です。 非圧縮性 Navier-Stokes 方程式の数値解法１：導入編 非圧縮性 Navier-Stokes 方程式の数値解法２：拡散項の陰解法（この記事） 非圧縮性 Navier-Stokes 方程式の数値解法３：陰解法の解き方比較 非圧縮性 Navier-Stokes 方程式の数値解法４：ソルバー精度検証 この記事の Livescript 版（MATLAB）は GitHub: 2D-Lid-Driven-Cavity-Flow, Part 21 に置いてありますので、実行しながら読みたい方はどうぞ。 この記事のゴール 導入編ではまずキャビティ流れを解いてみることに注力して、かなり荒い離散化を行いました。結果として、時間ステップサイズによっても変わりますが、ちょうどよい Reynolds 数（Re = 500程
非圧縮性 Navier-Stokes 方程式の数値解法１：導入編 - Qiita
はじめに Navier-Stokes 方程式を数値的に解きます。意外と簡単に流体シミュレーションできるんだな、、と思ってもらえると嬉しいです。間違い・分かりにくい箇所などはコメントください。随時更新していきます。 この記事の Livescript 版（MATLAB）はこちら GitHub: 2D-Lid-Driven-Cavity-Flow, Part11 に置いてありますので、実行しながら読みたい方はどうぞ。 この記事のゴール 「導入編」としてのゴールは 以下の図を作ること 部分段階法（Fractional Step Method）法を理解すること 離散化に際して Perot (1993) の境界条件の議論はあまり日本語で見かけないので、新しい情報になればと思います。 予告 途中で心変わりすると思いますがひとまず以下のような展開を予定しています。 非圧縮性 Navier-Stokes 方
スピン系のグラフ表現とImproved Estimator - Qiita
はじめに スピン系をモンテカルロシミュレーションする時、Swendsen-Wang法やWolff法など、Fortuin-Kasteleyn表現を用いたクラスターアップデートアルゴリズムを使うと非常に効率良くサンプリングができます。さらに、Fortuin-Kasteleyn表現が得られている時、普通にサンプリングするよりも分散が小さくなるImproved Estimatorを構成できます。 本稿では、Ising模型の磁化やビンダー比のImproved Estimatorの導出と、どのくらい分散が小さくなるのかの説明をします。 コードは https://github.com/kaityo256/mc/tree/master/improved_estimator に置いておきます。 Improved Estimatorとは サイコロとコインの例 Improved Estimatorを一言で言うと
数値計算の研究をしている学生が"数値計算に潜むとんでもないリスク"について話してみる - Qiita
筆者は「精度保証付き数値計算」という分野で研究をしている大学院生です． 「数値計算は分かるけど」「精度保証付き数値計算？ナニソレ？」という方がほとんどだと思います． 「精度保証付き数値計算」の研究自体は30年ほど前から盛んに行われていますが，世間に浸透しているとは言えない状況です． 自分の研究分野が世間に知られていないのは何か少し寂しい感じがするので「精度保証付き数値計算」を少しでも広めるべく記事を投稿することにしました．(シリーズ化するかも知れません) 本日は「精度保証付き数値計算」というワードだけでも覚えていただければ幸いです． 今回は"数値計算に潜むとんでもないリスク"に関してカジュアルにお話します． そして筆者の研究分野である「精度保証付き数値計算」の必要性を知ってもらえればなと思います． この記事を読み終える頃には計算機を信頼できなくなっているかも知れません(笑) ※不安を煽るこ
自動微分と区間演算で非線形関数の根を全て探索するアルゴリズム - Qiita
自動微分について調べていたら『精度保証付き数値計算の基礎』という本に出会いました。その中に非常に面白い応用が載っていたので紹介したいと思います。 allsol :: (RealFloat a, Ord a) => (forall b. Floating b => b -> b) -- 根を求める非線形関数 -> [Interval a] -- 探索する区間 -> [Interval a] -- 根が含まれている区間 これが今回紹介する、与えられた非線形関数$f$の与えられた区間における根、すなわち$f(x)=0$を満たすような$x$をただ一つ含む区間を全て探索してくれる関数（の型）です。実装を見る前にその威力を実際に使って確かめてみましょう。 > f x = (x - 1) * (x - 2) * (x - 3) > allsol f [-1e6...1e6] [ 0.9837477320
40行でPoisson方程式を解く - Qiita
こんにちは、やきつかです。 Haskellで配列を使ってみたいなーと思っていて、Poisson方程式を数値計算してみたら40行だったので1人で感動しちゃったので残しておきます。 Haskellの表現力の強さを知ってもらえればうれしいです。 注意 とりあえず実装したかったので、処理速度とかは一切考慮してないです。 同じ理由で、精度も可視化していい感じになることを目指しています。 Haskell、数値計算とも初学者なので、間違いとかよりよい作法がありましたら教えていただけると大変うれしいです。 はじめに この記事では、Poisson方程式の概要→差分方程式→アルゴリズム→Haskellで書いたコードという順で説明していきます。 Poisson方程式とは Poisson方程式は二階の偏微分方程式で、次の様に表されるもののことを言います。 $$\Delta\phi=f$$ ここで、$\Delta$
コンピュータ以前の数値計算(1) 三角関数表小史 -
現代の三角関数計算 三角関数の値を計算する方法として、現代人が素朴に思いつくのは (1)いくつかの角度に於ける値を事前に計算しておき、一般の場合は、それを補間した値を使う (2)Taylor展開の有限項近似 の二つの方法だと思う。Taylor展開を使う場合、角度をラジアン単位に変換する必要があるので、円周率を、ある程度の精度で知っていないといけない。 コンピュータ用に、もう少し凝ったアルゴリズムが使われることもある/あったらしいけど、今のコンピュータでは、(2)の方法が使われることが多い。例えば、Android(で採用されているBionic libc)では、アーキテクチャ独立な実装は、単純なTaylor展開を利用するものになっている。 https://android.googlesource.com/platform/bionic/+/refs/heads/master/libm/upst
格子計算プログラム生成言語Formuraを使ってみる その1 - Qiita
はじめに プログラマは計算機に対して、「別に凝ったをやれとは言わないが、自明なことはやってくれよ」と思うものだ。例えば規則格子における差分法、特に陽解法は、式と差分方式さえ決まれば、(最適化とか考えなければ)後は自明なコーディングになる。さらに、領域分割による分散メモリ並列を考えた時、(効率とか考えなければ)やはり自明なだけで面倒なコーディングをする必要がある。 こういう「自明かつ面倒なコーディング」をやらなければいけなくなった時、「自動化したい」と思うのは自然な発想だ。Formuraは、まさにそのような思想で作られた、格子計算プログラム生成言語(というかフレームワーク)である。 本稿では、このFormuraを使ってみる。 その1 インストールとコンパイルまで ←イマココ その2 一次元熱伝導方程式 その3 二次元熱伝導方程式 その4 反応拡散方程式(Gray-Scott系) インストール
閉殻Hartree-Fock法によるHeH+のエネルギー計算 - kivantium活動日記
量子化学計算を勉強するために新しい量子化学―電子構造の理論入門〈上〉を読んでいます。 サンプルコードがFortranだったので勉強がてらC言語に移植しました。 新しい量子化学―電子構造の理論入門〈上〉 作者:Attila Szabo,Neil S. Ostlund出版社/メーカー: 東京大学出版会発売日: 1987/07/01メディア: 単行本 理論 本で200ページくらいかけて説明されている内容をブログに書くのは大変なので省略します。 ちなみに、以下のツイートの「100ページかけて説明してる別の本」はこの本を指します。 難しい内容を易しく簡潔に説明しようとするのは大変で、難しい内容は難しいまま長い時間をかけて説明する必要があるのだなぁという学びがありました。 ある本で30ページで説明されてて何も分からなかった内容を100ページで説明してる別の本読んだら完全に理解して何が分からなかったのか
eom/bookの日本語版を作り始めた - Qiita
これは数値計算Advent Calendar23日目の記事です！クリスマスまでもう少しですね。 私は以前の記事 硬い微分方程式を解く でも紹介したeomというRustのライブラリを作っています。このライブラリは 常微分方程式、偏微分方程式（スペクトル法）の計算 力学系解析（リアプノフ指数、共変リアプノフベクトル）の計算 のためのライブラリです。 このライブラリは私の作ってる他のRustのライブラリに比べて数学的な側面が強いので、cargo docで生成されるライブラリとしてのリファレンスに加えて、数学的な側面を補足するbookを作っています。今回はAdCのためにこの日本語版を新たに書き始めました（できたとは言っていない）。 eom book/ja 数値計算の内容はbookを見てもらうとして、この記事ではせっかくなのでbookの作り方と管理方法についてまとめます。 mdbook 先程からbo
Vicsek model：二次元自己駆動粒子系と相転移 - Qiita
2018/12/21, 18:00 @goropikari さんにご指摘いただいた 粒子間距離の計算で周期境界を考慮できていない 近傍粒子の角度の計算で正負の打ち消しが生じている 問題を修正し，動画及び図を差し替えました． @goropikari さん，ありがとうございます． Julia Advent Calendar 2018，12/21日を担当します @Dolphin7473 です． 1. アクティブマターとVicsek model Julia advent calendarはみなさんJulia言語自体に関する興味深い記事を書いておられますが，「数値計算にJulia使ってみた」程度の記事が一つくらいあってもいいだろう，と思い執筆してみました．1 アクティブマターとは，ごく簡単に言うと 自走する粒子（群） です．バクテリアなどの微生物のような，各々が勝手に動きまわる物体（群）の総称で，こ
スペクトル要素法事始：Chebyshev多項式を用いて一次元Helmholtz方程式を解く - Qiita
スペクトル要素法：有限要素法とスペクトル法を通じた理解 TwitterではRunge-Kuttaに続いて有限要素法が流行っていますが，2018年の数値計算Advent calendarにもいくつか有限要素法を扱った記事がありました． @nqomorita さんの記事「有限要素法と反復法線形ソルバのはじめの一歩 [前編]」では，熱伝導方程式を題材に有限要素法の離散化がとてもわかりやすく解説されています． @t_kemmochi さんの記事「FreeFem++の紹介」はとてつもない大作で，有限要素法を用いた偏微分方程式の解析ソフトであるFeeeFem++のきわめて充実したチュートリアルとmatplotlibを用いた可視化，さらには付録として関数解析を用いた有限要素法の理論的な解説までがまとめられています． @muripo_life さんの記事「Galerkin法の数理についてまとめてみた」では
ルンゲ-クッタ法の安定性解析 - Qiita
はじめに この記事は数値計算 Advent Calendar 2018の１日目の記事です. ルンゲ-クッタ法によって求めた近似解の安定性を調べたいと思います. オイラー法 (1次のルンゲ-クッタ法) まずは1次のルンゲ-クッタ法である前進オイラー法(Forward Euler, FE)と後退オイラー法(Backward Euler, BE)の安定性を考えてみたいと思います. 常微分方程式の初期値問題
リリース、障害情報などのサービスのお知らせ
最新の人気エントリーの配信
j次のブックマーク
k前のブックマーク
lあとで読む
eコメント一覧を開く
oページを開く
GitHub - cantaro86/Financial-Models-Numerical-Methods: Collection of notebooks about quantitative finance, with interactive python code.
One Dimensional Heat Equation
GitHub - maplant/al-jabr: Extremely generic linear algebra libary for Rust
GitHub - cometscome/Julianotes: Julia language notes