はじめに

本書は、プログラミング言語Rustを用いた計算物理学の実践的な入門書です。数値計算の基礎から物理シミュレーションの実装までを体系的に解説し、Rustを用いた科学技術計算の基礎を身につけることを目的とします。

本書の目的

計算物理学は、現代の物理学研究において不可欠な手法です。解析的に解くことが困難な問題に対し、数値計算を用いることで具体的な解を得ることができます。一方、Rustは安全性とパフォーマンスを両立したプログラミング言語として注目を集めています。

本書は、これら2つの分野を結びつけ、Rustを用いて計算物理学の基本的な手法を実装することを通して、その原理の理解を深めることを目的としています。本書が、読者の皆様ご自身による物理シミュレーションの設計・実装、そして結果の分析の一助となることを目指します。

本書では、AI agent を用いた開発を推奨します。定型的な実装、入出力、テスト雛形、プロット用スクリプト、エラー調査、差分確認などは agent に任せることができます。一方で、生成されたコードの物理的・数値的な正しさは自動的には保証されません。本書では、Rust、数値計算、物理の基本概念をもとに、agent の出力を読み、実行し、検証する力を重視します。

想定読者

本書は、以下のような読者を想定しています。

• Rustを学びながら計算物理を実装したい方：変数、関数、所有権といったRustの基本概念を、必要に応じて公式ドキュメントを参照しながら学べる方を対象とします。高度なRustの知識は前提としませんが、簡単なプログラムを書いた経験があることを想定しています。

• 理工系の学部上級生から大学院修士初年程度の方：力学、電磁気学、熱力学・統計力学、量子力学といった基礎的な物理学の知識を前提とします。ただし、各章では必要に応じて関連する物理概念の簡単な復習を行います。

プログラミングや物理学の専門家である必要はありません。これらの分野を学びながら実践的なスキルを習得したい学生や研究者の皆様に最適な内容です。

Note

Rustを未習得の方へ

Rustに初めて触れる方は、本書と並行して、公式ドキュメントを必要に応じて参照してください。公式ドキュメントは非常に丁寧に記述されており、Rustの基本概念を体系的に学ぶことができます。

• The Rust Programming Language（英語版）
• The Rust Programming Language（日本語版）

特に、変数と可変性、所有権、構造体、列挙型、パターンマッチ、トレイトなどの章は、本書を読み進める上で重要な基礎となります。最初にすべてを完全に理解する必要はありません。演習で必要になった箇所を戻って確認する進め方で十分です。

本書の特徴

本書には、以下のような特徴があります。

• 実践重視のアプローチ：理論的な解説も行いますが、それ以上に、実際にコードを記述して実行し、結果を検証するプロセスを重視しています。各章で具体的な実装例を提示し、読者の皆様がご自身の手で試すことができるように構成しています。

• 体系的・網羅的な内容：数値計算の基礎に始まり、常微分方程式、偏微分方程式、モンテカルロ法、並列計算、可視化に至るまで、計算物理学の主要なトピックを幅広く網羅しています。

本書で扱う内容

本書は、大別して以下の4部構成です。

• 第1部 基礎編：Rustと計算物理学の関係、開発環境、計算機の基本モデル、 Rustで数値計算を書くための最小セット、多次元データ、ndarray、データ保存について解説します。

• 第2部 数値計算手法：科学計算の核となる技術を扱います。数値微分・積分、非線形方程式、線形代数、フーリエ解析、常微分方程式、偏微分方程式、モンテカルロ法などの数値解法を実装し、その特性の理解を深めます。

• 第3部 物理シミュレーション：実際の物理現象をシミュレートします。古典力学、流体力学、統計力学（イジング模型など）、量子力学といった、多様な物理分野における計算手法を解説します。

• 第4部 高度なトピック：より実用的なシミュレーション技術を扱います。並列計算による高速化、SIMD最適化、パフォーマンス測定などを通じて、実践的な技術を習得します。

本書の使い方

本書は、基本的に最初の章から順に読み進めることを推奨します。ただし、特定のトピックにご興味がある場合は、必要な章から読み始めることも可能です。第1部（基礎編）は、それ以降の章における前提知識となる項目が多いため、一度目を通されることを推奨します。

各章は、概ね以下の構成で展開されます。

1. 概念の理解：トピックの物理的・数学的背景を簡潔に解説します。
2. 実装：Rustによる具体的な実装方法を示します。
3. 結果保存と可視化：コードを実行し、得られた結果を保存して確認します。
4. 考察：結果の解釈や、手法の特性について考察します。

掲載されているコードは、ご自身の手で実行してみることを強く推奨します。エラーへの対処や、期待した結果が得られない場合の試行錯誤は、学習における重要なプロセスです。こうした経験を通じて、より深い理解を得ることができます。

本書について

本書は、櫻庭が卒業研究として作成したRust計算物理学チュートリアルをもとに、 品岡寛が全体構成、科学的内容、演習、AI agent を用いた開発方針を拡張したものです。

執筆にあたり、AI（大規模言語モデル）を活用していますが、生成された内容はすべて著者によって事実確認およびコードの動作検証が行われています。AIは執筆を効率化するための補助的なツールとして用いており、内容の正確性については著者が責任を負います。

本書が、これから計算物理学を学ばれる皆様や、Rustを用いた科学技術計算にご興味をお持ちの方々にとって、有益な一助となれば幸いです。

フィードバックと貢献

本書の内容には万全を期しておりますが、誤字脱字、理論的な誤り、 あるいはコードの不具合などが含まれている可能性がございます。 もし何かお気付きの点や改善のご提案がございましたら、 GitHubリポジトリ へのIssueの作成、またはPull Requestを送付いただけますと幸いです。 ご報告は日本語で問題ございません。 皆様からのフィードバックは、本書の品質向上に不可欠です。

また、「この説明が理解しにくい」「このようなトピックも扱ってほしい」といったご意見も歓迎いたします。読者の皆様と共に、より良い教材を構築していければ幸いです。

ライセンス

本書は、CC0 1.0 Universal (CC0 1.0) Public Domain Dedicationライセンスの下で公開されています。著作者は、著作権法上認められる最大限の範囲で、本書に対する著作権および関連する権利を放棄しています。本書の内容は、教育、研究、商用を問わず、自由に利用、改変、再配布が可能です。

それでは、Rustによる計算物理学の世界をお楽しみください。

第1章 Rustと計算物理学

本章は、プログラミング言語Rustを用いた計算物理学の学習を開始するための導入部です。 本章では、なぜ本書でRustを使うのか、どのように開発環境を整えるのか、計算結果をどのように保存・可視化するのかを確認します。

本書では、Rustを「人間がすべて手で書くための言語」としてだけではなく、AI agent が生成・修正したコードをコンパイラ、テスト、差分確認の中で検査しやすい言語として扱います。重要なのは、コードを書く速度そのものではなく、物理モデル、数値計算法、結果保存、検証を再現可能な形で進めることです。

本章の構成

本章は、以下の内容で構成されます。

• なぜRustなのか？ Rustの型システム、所有権、Cargo、テストが、計算コードとAI agent の出力を検査する足場になることを説明します。

• 開発環境のセットアップ Rust toolchain、cargo、gnuplot、エディタ、推奨するAI coding agent の位置づけについて確認します。

• 数値計算結果の可視化 Rustで計算し、結果をCSVやバイナリ形式で保存し、保存済みデータから図を作る方針を説明します。表形式データではgnuplotを標準の可視化手段として使います。

• 本書の使い方 各章を読む前に確認すべき問題設定、入力と出力、検証方法、AI agent に任せる作業と人間が確認する作業を整理します。

本章を終えることで、Rustによる数値計算を始めるための環境と、本書を読み進めるための基本的な作業方針が整います。

なぜRustを用いるのか

従来、計算物理学の分野では、Fortran、C/C++、Python、Julia、MATLABといったプログラミング言語が広く利用されてきました。これらの言語にはそれぞれ長所があり、特にPythonは科学技術計算におけるエコシステムが非常に充実しています。

AI agent を用いてコードの作成や修正を補助する時代には、言語選択の基準が少し変わります。人間がすべてのコードを手で書きやすいかだけでなく、人間が読みやすいか、生成されたコードを検査しやすいか、誤った変更をコンパイラやテストで早く止められるかが重要になります。

本書でRustを採用する理由は、AI agent がRustコードを常に正しく生成できるからではありません。むしろ、Rustの型システム、所有権、借用、明示的な可変性、Cargo、cargo testが、生成されたコードを検査するための強い足場になるからです。

AI agent 時代におけるRustの利点

1. 生成コードをコンパイラで制約できる

計算物理のコードでは、配列、状態変数、乱数生成器、作業用バッファ、出力データなど、多くのデータを扱います。これらのデータをいつ変更してよいのか、どこで共有してよいのかが曖昧になると、非常に追跡しにくいバグが生じます。

C/C++では、aliasing、メモリの寿命、可変な共有、データ競合に起因する問題が、実行時のまれな条件で初めて現れることがあります。Rustでは、所有権、借用、ライフタイム、Send / Sync などの仕組みにより、この種の問題の多くをコンパイル時に検査できます。

これは、AI agent が生成したコードを無制限に信用するのではなく、コンパイラとテストによるハーネスの中で扱うという意味で重要です。cargo checkで型や借用の問題を早く見つけ、cargo testで小さい計算部品を検証し、git diffで変更内容を確認する、という流れを作りやすくなります。

ただし、コンパイルが通ることは、物理的・数値的に正しいことを意味しません。数式、単位、境界条件、刻み幅、乱数 seed、保存量、収束性などは、読者自身が問題設定に立ち返って検証する必要があります。

2. 見た目の複雑さはレビュー不能性ではない

Rustのコードは、PythonやJuliaのコードに比べて複雑に見えることがあります。型、参照、mut、Result、エラー処理、モジュール境界など、数値計算の式そのもの以外の記述が多いためです。

しかし、その多くはコンパイラに渡す制約や、コードの意図を明示するための情報です。ロジックを追うときは、まず物理量の更新、ループ、境界条件、保存する値に注目し、型注釈やエラー処理の定型部分を一旦脇に置くことができます。むしろ、どの値が可変か、どこで所有権が移るか、どの失敗を扱うかが明示されているため、曖昧なコードよりもレビューしやすい場面があります。

AI agent を使う場合は、レビューの補助も依頼できます。例えば、「この関数を疑似コードにして」「物理量の更新ロジックだけを抜き出して」「この Result はどの失敗を表しているか説明して」と指示すれば、Rust特有の記述とアルゴリズム本体を分けて確認できます。

3. 学習コストの意味が変わる

Rustは、学習コストが高い言語として知られてきました。所有権、借用、ライフタイム、トレイトなどは、他の多くの手続き型言語にはない概念です。従来は、細かい記法やコンパイラエラーへの対応そのものが、Rustを使う上で大きな負担になりがちでした。

しかし、AI agent を用いる開発では、コードを最初から最後まで人間が手で書くことだけが中心ではありません。定型的な構造、入出力、テストの雛形、エラー処理、モジュール分割などは、AI agent に任せることができます。その分、人間は、物理モデル、数値計算法、境界条件、検証方法、結果の解釈に集中しやすくなります。

このとき、Rustの厳しさは単なる欠点ではなくなります。コンパイラが細かい誤りを指摘し、AI agent が修正案を出し、人間が物理・数値計算として妥当かを判断する、という分担が作りやすいからです。

4. エコシステム不足の見方も変わる

Pythonの科学計算エコシステムは、NumPy、SciPy、Matplotlib、pandasなどを中心に非常に成熟しています。一方、Rustの科学技術計算エコシステムは発展途上であり、特定の数値計算アルゴリズムや専門分野に特化したクレートが存在しない場合もあります。

これは現実的な課題です。ただし、AI agent を用いると、薄い wrapper、入出力コード、テスト、保存済みデータからのプロットなどは、以前よりも短い時間で用意しやすくなります。すべての機能が既にクレートとして揃っているからRustを使うのではなく、不足している部分を検査可能な形で補いやすいことが重要になります。

従来の言語との関係

Rustは、既存の言語をすべて置き換えるものではありません。目的に応じて、以下のような違いがあります。

• Fortran：数値計算やHPCの分野で長い実績があり、既存の数値ライブラリ資産に強みがあります。一方で、テンソルネットワークのように複雑なデータ構造、抽象化、テスト分割、依存関係管理を伴うアルゴリズムでは、設計や保守の負担が大きくなりやすい面があります。
• C/C++：高速な実行が可能で、既存資産も豊富です。一方で、メモリ管理、aliasing、未定義動作、データ競合などの問題をプログラマが慎重に管理する必要があります。
• Python：探索的な解析、可視化、既存ライブラリの利用に非常に強い言語です。一方、言語単体の実行速度や静的検査には限界があり、高速化が必要な箇所はC、Fortran、Rustなどで実装された外部ライブラリに依存することが多くなります。
• Julia：高水準な記述と高い実行性能を両立しようとする言語で、数値計算のプロトタイピングにも適しています。一方、Rustのように所有権や借用を用いてメモリ共有やデータ競合をコンパイル時に厳密に制約する設計ではありません。

本書では、Rustを「最も簡単に書ける言語」としてではなく、計算コードを安全に育て、AI agent の出力を検査し、再現可能な形で残しやすい言語として扱います。

Rustの具体的な強み

安全性とパフォーマンスの両立

Rustは、所有権システムによりコンパイル時にメモリ安全性を保証しつつ、C/C++に近い実行性能を目指せる言語です。ガベージコレクタに依存しないため、実行中の予期しにくい停止を避けやすく、長時間実行する計算や性能が重要な計算に向いています。

データ競合のない並列処理

現代の計算物理学では、マルチコアCPUやGPUといった並列計算資源の活用が重要です。Rustの型システムは、スレッド間で共有してよいデータと共有してはいけないデータをコンパイル時に検査します。

例えば、データ並列ライブラリであるRayonを利用すると、イテレータを並列化しながら、データ競合を起こしにくい形でコードを書くことができます。

Cargoによる再現可能な開発環境

Rustでは、cargo newによるプロジェクト作成、cargo addによる依存関係の追加、cargo buildによるビルド、cargo testによるテスト実行といった一連の作業が、cargoという単一のツールで完結します。

科学技術計算では、どの依存関係で、どの入力条件で、どのコードを実行したのかを後から確認できることが重要です。Cargo.tomlとCargo.lockにより依存関係を管理し、テストと実行コマンドを残すことで、計算の再現性を高められます。

他言語との連携

Rustからは、Fortran、C、C++で実装された既存の数値計算ライブラリを呼び出すこともできます。過去の資産をすべて捨てるのではなく、必要に応じて利用できる点も実用上の利点です。

使い分け

Rustにも向き不向きはあります。学習コストは依然として存在しますし、Pythonほど科学計算ライブラリが網羅的に揃っているわけではありません。

本書では、Rustで計算を行い、表形式データはCSV、多次元配列や大きなデータはHDF5やNumPy形式、実行条件はmetadataとして保存し、プロットは保存済みデータから作成する流れを基本にします。この方法では、入力パラメータ、乱数 seed、実行条件、出力ファイルを明示的に残せるため、結果の再現性や検証がしやすくなります。入出力やプロット用の定型的なコードは、AI agent に自然言語で指示して任せることができます。

迅速なデータ分析、使い捨てのスクリプト、既存Pythonライブラリに強く依存した解析では、PythonやJuliaを使う方が自然な場合もあります。一方でRustは、長時間実行する計算、性能やメモリ効率が重要な計算、並列化したい計算、長期的に保守したい計算コードで特に有用です。

まとめ

Rustは、安全性、パフォーマンス、並列処理、再現可能な開発環境という点で、計算物理学にとって魅力的な選択肢です。AI agent 時代には、さらに、生成コードをコンパイラとテストのハーネスの中で扱えることが重要になります。

AI agent は、コードを速く書く助けにはなりますが、物理的・数値的な正しさを自動的に保証するものではありません。本書では、Rustの厳密な仕組みを利用しながら、読者自身が問題設定、実装、検証、結果の解釈を行えるようになることを目指します。

開発環境の構築

本章では、Rustを用いて計算物理学の研究・開発を開始するための開発環境を構築します。本書は、Rustの基本文法を既に習得されている読者を対象としているため、多くの方はRustをインストール済みであると想定されます。該当する方は、本節を読み飛ばしていただいて問題ありません。

Rustをまだインストールされていない方や、環境を再確認したい方のために、標準的な構築手順を以下に解説します。

Rustのインストール

Rustをまだインストールされていない方は、Rust公式サイトの案内に従ってインストールを実行してください。公式サイトでは、お使いのOS（Windows, macOS, Linux）に応じた最新のインストール手順が提供されています。

Rust公式のツールチェーン管理ツールであるrustupを使用することで、Rustコンパイラのバージョン管理や、関連ツールの更新を容易に行うことができます。

インストールの確認

Rustが正しくインストールされているかを確認するため、ターミナルで以下のコマンドを実行してください。

rustc --version
cargo --version

これらのコマンドがそれぞれのバージョン情報を表示すれば、インストールは成功です。

• rustc: Rustのコンパイラ
• cargo: Rustのビルドツール兼パッケージマネージャ

本書では、プロジェクトの作成や管理にcargoを主に利用します。

Rustのバージョンについて

本書では、Rustの安定版（stable）を使用することを推奨します。Rustは後方互換性が重視されているため、多くの場合、少し古いバージョンでも本書のコードは問題なく動作します。

古いバージョンのRustを使用していてコンパイルエラーが発生した場合は、以下のコマンドで更新してください。

rustup update

特定のバージョンを使う必要がある場合は、rust-toolchain.tomlなどでプロジェクトごとに固定できます。ただし、本書を読み進める上では、まずstableを使えば十分です。

可視化ツールの確認

本書では、表形式のデータをCSVに保存し、gnuplotで可視化する流れを標準とします。gnuplotが利用できるかは、次のコマンドで確認できます。

gnuplot --version

インストールされていない場合は、お使いのOSのパッケージマネージャや公式サイトの案内に従って導入してください。

Pythonのmatplotlib、HDF5、NumPy形式（.npy、.npz）は必須ではありません。複雑なデータ加工、多次元配列、大きなデータを扱う章や発展課題で、必要に応じて利用します。

開発環境の動作確認

次に、簡単なプログラムを作成・実行し、開発環境が正しく動作することを確認します。ここでは、本書のテーマに合わせて、簡単な数値計算を行うプログラムを作成します。

プロジェクトの作成

まず、新しいプロジェクトを作成します。任意のディレクトリで以下のコマンドを実行してください（physics_testの部分は任意のプロジェクト名に変更可能です）。

cargo new physics_test
cd physics_test

cargo newは、新しいRustプロジェクトの雛形を生成するコマンドです。physics_testという名前のディレクトリが作成され、その中に以下のようなファイル構造が生成されます。

physics_test/
├── Cargo.toml     # プロジェクトの設定ファイル（依存関係など）
└── src/
    └── main.rs    # メインのソースファイル

簡単な数値計算プログラムの作成

src/main.rsをエディタで開き、以下の内容に書き換えてください。このプログラムは、自由落下する物体の位置と速度を時系列で計算するものです。

fn main() {
    // 物理定数
    let g = 9.8; // 重力加速度 [m/s^2]
    let v0 = 0.0; // 初速度 [m/s]
    let h0 = 100.0; // 初期高度 [m]

    println!("自由落下シミュレーション");
    println!("初期高度: {} m", h0);
    println!("重力加速度: {} m/s^2", g);
    println!();

    // 時刻ごとの位置と速度を計算
    println!("時刻[s]  高度[m]  速度[m/s]");
    println!("================================");

    for i in 0..=10 {
        let t = i as f64 * 0.5; // 0.5秒刻み
        let v = v0 - g * t; // 速度 v = v0 - gt
        let h = h0 + v0 * t - 0.5 * g * t * t; // 位置 h = h0 + v0*t - (1/2)*g*t^2

        if h >= 0.0 {
            println!("{:6.1}  {:7.2}  {:9.2}", t, h, v);
        } else {
            // 地面に到達した場合、正確な到達時刻と速度を算出して終了
            let t_impact = (v0 + (v0 * v0 + 2.0 * g * h0).sqrt()) / g;
            let v_impact = v0 - g * t_impact;
            println!("{:6.2}  {:7.2}  {:9.2} (地面到達)", t_impact, 0.0, v_impact);
            break;
        }
    }
}

このプログラムは、以下の物理法則に基づき実装されています。

• 速度：
• 位置：

ここで、 は重力加速度、 は初速度、 は初期高度です。

プログラムの実行

プロジェクトのルートディレクトリ（Cargo.tomlが配置されているディレクトリ）で、以下のコマンドを実行してください。

cargo run

初回の実行ではプロジェクトのコンパイルが行われるため、多少時間を要します。コンパイルが完了するとプログラムが自動的に実行され、以下のような出力が表示されます。

自由落下シミュレーション
初期高度: 100 m
重力加速度: 9.8 m/s^2

時刻[s]  高度[m]  速度[m/s]
================================
   0.0   100.00       0.00
   0.5    98.78      -4.90
   1.0    95.10      -9.80
   1.5    88.97     -14.70
   2.0    80.40     -19.60
   2.5    69.38     -24.50
   3.0    55.90     -29.40
   3.5    39.97     -34.30
   4.0    21.60     -39.20
   4.5     0.77     -44.10
  4.52     0.00     -44.27 (地面到達)

このような出力が得られれば、開発環境は正しく構築されています。

コードの概要

上記のコードでは、Rustの基本的な機能をいくつか利用しています。

• 変数の宣言: letキーワードで変数を宣言します。Rustの変数はデフォルトでイミュータブル（不変）です。
• 浮動小数点数: f64型（64ビット浮動小数点数）を利用しています。多くの場合、型はコンパイラによって推論されるため、明示的な記述は不要です。
• ループ: forループを用いて、0.5秒ごとに時間ステップを進めています。
• 型変換: i as f64という構文で、ループカウンタの整数iを浮動小数点数に変換しています。
• 条件分岐: if文を用いて、物体が地面に到達したか否かを判定しています。
• 書式付き出力: println!マクロと書式指定子（例: {:6.1}）を使い、整列された出力を生成しています。

このプログラムは非常に単純ですが、計算物理学の基本的な要素（物理法則の数式化、時間発展シミュレーション、結果の出力）を含んでいます。本書では、これらの要素をさらに発展させ、より複雑な物理現象のシミュレーションを扱います。

エディタとIDE

Rustのコードを記述するためのエディタや統合開発環境（IDE）は、ご自身の好みに応じて選択いただけます。多くの主要なエディタには、Rustをサポートする拡張機能が提供されています。

どのエディタを選択する場合でも、rust-analyzerのサポートを有効にすることを強く推奨します。rust-analyzerは、Rust公式の言語サーバープロトコル（LSP）実装であり、以下のような強力な開発支援機能を提供します。

• コード補完
• エラーチェックとハイライト
• 定義元へのジャンプ
• リファクタリング

これらの機能により、開発効率が大幅に向上します。

AI coding agent の利用（推奨）

本書の演習では、AI coding agentの利用を推奨します。例えば、Codex、OpenCode、Claude Code などのツールを用いることで、コードの作成や修正、入出力、テスト雛形、プロット用スクリプト、エラー原因の調査、差分の確認を任せることができます。

ただし、特定のツールを必須とはしません。各ツールのインストール方法、料金、利用条件は変わる可能性があるため、本書では詳細な手順や比較は扱いません。利用する場合は、それぞれの公式情報を確認してください。

実践的な導入資料として、CompPhysHack 2026 の講義資料と、 AtelierArith による coding agents のハンズオン資料も参考になります。 後者では、Claude Code や Cursor などの coding agent を用いた開発の進め方が、 演習を通して紹介されています。

• CompPhysHack 2026
• CompPhysHack2026HandsOn by AtelierArith

AI coding agent は、演習用プロジェクトや本書のリポジトリのルートディレクトリで起動し、 Cargo.toml、srcディレクトリ、テスト、差分を参照できる状態にします。 agent に任せた変更は、cargo checkやcargo testで確認し、 必要に応じてgit status、git diff、git diff HEADで変更内容を確認してください。

Rustのコンパイラとテストは、AI coding agent が生成したコードを検査するためのハーネスとして機能します。一方で、コンパイルが通ることは、物理モデル、数値計算法、単位、境界条件、乱数 seed などが正しいことを意味しません。これらは、各章で扱う理論と検証方法に立ち返って確認する必要があります。

次のステップ

以上で、開発環境の構築と、簡単な物理シミュレーションの実行が完了しました。しかし、画面への出力だけでは、結果を直感的に理解することは困難です。

次節では、シミュレーション結果をグラフとして可視化する手法について解説します。

トラブルシューティング

もしインストールや実行時に問題が発生した場合は、付録C: デバッグとトラブルシューティングをご参照ください。そこでは、よくある問題とその解決策が記載されています。

また、Rust公式のドキュメントも非常に充実しており、問題解決の一助となります。

• The Rust Programming Language（日本語版）：Rustの公式入門書（通称 The Book）
• Rust By Example（日本語版）：豊富な実例を通じてRustを学ぶ
• The Cargo Book：Cargoの公式ドキュメント

数値計算結果の可視化

数値計算では、計算結果を図として確認することが重要です。ただし、本書の中心は可視化ライブラリの使い方ではなく、物理モデル、数値計算法、Rustによる実装、そして結果の検証です。

そのため本書では、まず次の流れを標準とします。

1. Rustで計算する。
2. 結果をCSVやバイナリ形式のファイルに保存する。
3. 表形式データはgnuplotで可視化する。
4. 図を見て、解析解、保存量、収束性、境界条件などを確認する。

gnuplotはテキストファイルでプロット手順を管理できるため、表形式データの確認には十分です。多次元配列では、必要な断面、射影、集計量を取り出して可視化します。複雑なデータ加工や論文用の細かい図が必要になった場合だけ、Pythonのmatplotlibなどを使います。

基本方針

本書では、計算と可視化を分離します。Rustのプログラムは、画面への出力だけでなく、後から確認できる形でデータを保存します。

この方針には、次の利点があります。

• 入力パラメータ、乱数seed、実行条件、出力ファイルを後から確認できる。
• 同じデータから何度でも図を作り直せる。
• 図の作成手順をplot.gpのような短いテキストファイルとして管理できる。
• AI agentに、入出力コードやプロット用スクリプトの作成を任せやすい。

可視化は、計算結果を検査するための手段です。図がきれいに描けることと、数値計算が正しいことは同じではありません。

CSVファイルへの出力

CSVは多くのツールで読み込めるため、時系列や表形式の小さなデータに適しています。次の例では、時刻t、位置x、速度vをoutput.csvに保存します。

use std::fs::File;
use std::io::Write;

fn main() -> std::io::Result<()> {
    let mut file = File::create("output.csv")?;

    writeln!(file, "t,x,v")?;

    for i in 0..=100 {
        let t = i as f64 * 0.1;
        let x = t.sin();
        let v = t.cos();
        writeln!(file, "{t},{x},{v}")?;
    }

    Ok(())
}

実行すると、次のようなCSVファイルが生成されます。

t,x,v
0,0,1
0.1,0.09983341664682815,0.9950041652780258
0.2,0.19866933079506122,0.9800665778412416

実際のシミュレーションでは、必要に応じてmetadata.jsonのような別ファイルに、パラメータ、乱数seed、格子サイズ、刻み幅、実行日時なども保存するとよいでしょう。

多次元配列とバイナリ形式

格子上の場、画像、波動関数、相関関数などの多次元配列を扱う場合、CSVはあまり適していません。ファイルサイズが大きくなりやすく、形状やデータ型の情報も別に管理する必要があるためです。

このような場合は、HDF5やNumPy形式（.npy、.npz）などのバイナリ形式を検討します。これらは広く使われており、Python、Julia、C/C++、Fortranなど複数の言語から読み書きしやすい形式です。

したがって、データ形式は次のように使い分けます。

• 時系列や少数列の表形式データ: CSV
• 多次元配列や大きなデータ: HDF5、NumPy形式などのバイナリ形式
• 実行条件やパラメータ: JSON、TOML、YAMLなどのメタデータファイル

gnuplotによる可視化

gnuplotでは、プロット手順をテキストファイルとして保存できます。例えば、次の内容をplot.gpとして保存します。

set terminal svg size 900,600
set output "plot.svg"
set datafile separator ","
set key autotitle columnhead

set xlabel "t"
set ylabel "value"
set grid

plot "output.csv" using 1:2 with lines linewidth 2, \
     "" using 1:3 with lines linewidth 2

次のコマンドで図を生成します。

gnuplot plot.gp

この例では、output.csvの1列目を横軸、2列目と3列目を縦軸として描画し、plot.svgを出力します。CSVのヘッダー行は凡例として使われます。

gnuplotスクリプトは短く、差分も読みやすいため、agentic codingと相性がよいです。例えば、AI agentには次のように依頼できます。

output.csvには t, energy, error が入っています。
energyの時間変化とerrorの対数プロットを別々のSVGに保存するgnuplotスクリプトを書いてください。
軸ラベル、grid、出力ファイル名も入れてください。

ただし、生成された図を確認するのは人間の役割です。軸、単位、列の対応、期待される保存量や収束性が正しいかを確認してください。

matplotlibを使う場合

Pythonのmatplotlibは、データ加工や複雑な図を作る場合に有用です。例えば、複数ファイルの集計、フィッティング、統計量の計算、論文用の細かな体裁調整を行う場合にはmatplotlibを使うと便利です。

import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_csv("output.csv")

fig, ax = plt.subplots()
ax.plot(df["t"], df["x"], label="x")
ax.plot(df["t"], df["v"], label="v")
ax.set_xlabel("t")
ax.set_ylabel("value")
ax.grid(True)
ax.legend()
fig.savefig("plot.png", dpi=150)

一方で、単純な折れ線グラフや散布図であれば、まずgnuplotで十分です。可視化のためにPython環境を必須にしないことを、本書の標準方針とします。

Rust製ライブラリを使う場合

Rustだけで図の生成まで完結させたい場合は、plottersを使うこともできます。また、3D可視化やインタラクティブな表示が必要な場合には、three-dなどのグラフィックスライブラリが選択肢になります。

ただし、これらは本書の標準経路ではありません。Rust側に可視化処理を深く組み込むと、数値計算の本体と図の作成が混ざりやすくなります。まずは計算結果を保存し、外部ツールで可視化する構成を基本にしてください。

まとめ

本書では、次の方針を採ります。

Important

• Rustは計算と結果保存を担当する。
• 標準の可視化は、表形式データならCSVとgnuplotを使う。
• 多次元配列や大きなデータは、HDF5やNumPy形式などのバイナリ形式で保存する。
• matplotlibは、複雑なデータ加工や図が必要な場合の補助として使う。
• plottersやthree-dは、Rustだけで完結させたい場合や3D可視化が必要な場合の発展的な選択肢とする。
• 図は検証の入口であり、物理的・数値的な正しさは別途確認する。

次節では、本書全体の構成と学習の進め方について解説します。

本書の使い方

本節は、本書を最大限に活用し、計算物理学とRustプログラミングの技術を効率的に習得するための指針を示すものです。本書は、理論の学習、コードの実装、結果の保存、そして検証と考察というサイクルを通じて、実践的な能力を養うことを目的として設計されています。

基本的な進め方：通読

計算物理学や数値計算の全体像を体系的に学習したい方には、第1章から順に読み進めることを強く推奨します。前の章で解説した概念や実装が、後の章の基礎となっていることが多いため、この方法が最もスムーズに理解を深めることができます。

もちろん、興味のある章から読み進めても構いませんが、第1部（基礎編） は本書全体の共通の基礎となるため、一度目を通しておくことを推奨します。

実装に入る前に確認すること

本書の中核は、実際にコードを動かし、結果を検証する過程にあります。ただし、計算物理のコードでは、いきなり実装から始めると、何を確認すべきかが曖昧になりがちです。

各章の演習では、まず次の点を確認してください。

1. 問題設定：何を計算したいのか。
2. 入力と出力：パラメータ、初期条件、保存する結果は何か。
3. 仮定と境界条件：どの近似や境界条件を使うのか。
4. 数値計算法：どのアルゴリズムを使い、どの量を更新するのか。
5. データ構造：配列、状態変数、乱数生成器、結果保存をどう表すのか。
6. 検証方法：解析解、保存量、収束性、極限ケース、既知の結果とどう比較するのか。

AI agent を使う場合でも、この確認を省略しないでください。agent に実装を任せる場合も、人間が問題設定と検証方針を持っていることが重要です。

サンプルコードの実践ガイド

本書のサンプルコードは、以下の手順で実行することを想定しています。

1. 新しいプロジェクトを作成する ターミナルで cargo new コマンドを用いて、演習用のプロジェクトを作成します。

   cargo new physics-exercise
   cd physics-exercise
   

2. コードを配置する 任意のテキストエディタで src/main.rs や src/lib.rs を開き、本書に掲載されているコード例を入力します。

3. 依存関係を追加する Cargo.toml ファイルを開き、コードが必要とするクレート（ライブラリ）を [dependencies] セクションに追加します。バージョンは、本書の表記に合わせることを推奨します。

   [dependencies]
   csv = "1"
   # 章で必要なクレートだけを追加
   

4. コンパイルと実行 ターミナルで cargo run を実行してください。Cargoが自動的に依存クレートをダウンロードし、コードをコンパイルして実行します。

   cargo run
   

5. 確認する 実装後は、少なくとも次のコマンドを確認します。

   cargo test
   git status
   git diff
   git diff HEAD
   

   cargo test は、実装した関数や小さい検証問題が期待通りに動くかを確認するために使います。 git status は変更されたファイルの確認、git diff はまだstageしていない差分の確認、 git diff HEAD は最終commitからの変更全体の確認に使います。

AI agent を使う標準的な流れ

本書では、AI agent の利用を推奨します。ただし、いきなり「この章のコードを実装して」と頼むのではなく、まず問題設定と検証方針を整理させるとよいでしょう。

例えば、次のような順序で進めます。

1. note を作らせる モデル定義、入力、出力、仮定、境界条件、数値計算法、検証方法をまとめさせます。

2. implementation plan を作らせる 関数境界、module 構成、テスト、保存する結果、metadata、実行コマンドを含む計画を作らせます。

3. 実装を任せる 定型的なコード、入出力、テスト雛形、プロット用コード、モジュール分割は agent に任せることができます。ただし、物理モデルや数値計算法として妥当かどうかは、人間が確認します。

4. 検査する cargo test、保存された結果、プロット、git diff HEAD を確認します。コンパイルが通ることは、物理的・数値的に正しいことを意味しません。

5. 具体的にレビューさせる 「正しいですか」とだけ聞くのではなく、次のように確認対象を指定します。

   • 解析解と比較してください。
   • 保存量が時間発展でどの程度保たれているか確認してください。
   • 刻み幅を変えたときの収束性を確認してください。
   • 境界条件の実装を確認してください。
   • random seed、入力パラメータ、出力ファイル名が保存されているか確認してください。
   • この関数を疑似コードにして、物理量の更新ロジックだけを説明してください。

結果保存とプロット

本書では、Rustで計算を行い、表形式データはCSV、多次元配列や大きなデータはHDF5やNumPy形式、実行条件はmetadataとして保存し、プロットは保存済みデータから作成する流れを基本にします。

この方法には、次の利点があります。

• 入力パラメータ、乱数 seed、実行条件、出力ファイルを後から確認できる。
• 計算と可視化を分離できる。
• 同じ保存データから、何度でもプロットを作り直せる。
• AI agent に、保存処理やプロット用コードを自然言語で任せやすい。

単に画面に数値を表示するだけでなく、後から検証できる形で結果を残す習慣をつけてください。

コードを用いた実験の推奨

コードの正常な動作を確認した後は、パラメータや初期条件を変更し、「実験」を行うことを強く推奨します。これは、内容の理解を深める上で極めて効果的な学習方法です。

• パラメータの変更: シミュレーションのステップ数、物理定数、初期値などを変更し、結果がどのように変化するかを観察してください。
• 題材の応用: sin(x) をプロットする代わりに、別の関数（例: cos(x) やガウス関数）をプロットしてみてください。
• エラーからの学習: エラーは絶好の学習機会です。コンパイラが示すメッセージを注意深く読み、問題の原因を突き止めるプロセスを通じて、Rustへの理解が格段に深まります。

もしコードが期待通りに動作しない場合や、完成版のコードを確認したい場合は、本書の公式リポジトリをご参照ください。ただし、まずはご自身で原因を調べ、必要に応じてAI agent に差分やエラーの説明を求めることを推奨します。

各章の構成とヒント

本書の多くの章や節は、概ね以下の構成になっています。

1. 理論的背景: トピックの背後にある物理的・数学的な概念を解説します。
2. 実装: Rustによる具体的な実装方法を示します。
3. 実行と結果保存: コードを実行し、得られた結果を保存します。
4. 可視化と考察: 保存済みデータからグラフを作成し、アルゴリズムの特性や物理的な意味を考察します。
5. 検証: 解析解、保存量、収束性、境界条件、metadata などを確認します。

章の最後には、学んだ知識をさらに深めるための 「発展課題」 を設けていることがあります。これらには必ずしも唯一の正解があるわけではなく、読者の皆様が自由な探求を行うための指針として提示されています。ぜひ挑戦されることを推奨します。

本書の表記ルール

本書では、情報を明確に伝達するため、いくつかの特別な表記法を用いています。

Note

NOTE (補足)

補足情報や、関連する豆知識などを記載しています。

Important

IMPORTANT (重要)

特に重要な概念や、必ず留意していただきたい事柄を示しています。

Warning

WARNING (警告)

注意を怠るとエラーや予期せぬ結果に繋がりかねない、重要な注意点を示しています。

以上で、本書の学習を開始するための準備はすべて整いました。

次章から、Rustによる数値計算の学習を開始します。

計算機の基本モデル

計算物理のコードは、数式だけでなく、実際の計算機の上で動きます。 同じ計算量のアルゴリズムでも、メモリ上のデータの並びやアクセス順序によって、 実行時間が大きく変わることがあります。

本章では、Rustに入る前に、数値計算コードを読むための計算機モデルを確認します。 ここで扱う内容は Rust 固有ではありません。C/C++、Fortran、Python/NumPy でも重要です。

• CPU、メモリ、ストレージの役割
• latency と bandwidth
• cache と cache line
• 1次元データの連続アクセスと飛び飛びアクセス
• stack と heap

本章の目的は、CPU architecture を詳しく学ぶことではありません。 以降の章で Vec<f64>、多次元配列、ndarray、PDE、行列演算、 SIMD、並列化を扱うときに、なぜデータ配置が重要なのかを理解することです。

本章の構成

• CPU・メモリ・ストレージ 計算がどこで行われ、データがどこに置かれるかを整理します。

• latency・bandwidth・cache line メモリアクセスの遅さ、帯域、cache line の直感を確認します。

• 連続アクセスと1次元データ Vec<f64> のような1次元データで、連続アクセスが性能に効く理由を見ます。

• stackとheap 小さい値と大きな配列がどこに置かれるかを確認します。

本章を終えると、以降の章で出てくる配列操作や性能の話を、 プログラミング言語だけではなく計算機の動きとして読めるようになります。

CPU・メモリ・ストレージ

計算機を大まかに見ると、次の3つの場所があります。

• CPU: 命令を実行し、四則演算や分岐を処理する。
• メモリ: 実行中のプログラムが使うデータを置く。
• ストレージ: ファイルとしてデータを長期保存する（SSD、HDDなど）。

数値計算では、CPUがすべての時間を計算に使っているとは限りません。 大きな配列を扱う場合、CPUはメモリからデータが届くのを待っていることがあります。 このような場合、コードの速さは浮動小数点演算の回数だけでは決まりません。

大まかな接続

現代のPCやサーバのCPUには、通常、複数のcoreがあります。 各coreが命令を実行し、複数coreを使うと並列計算につながります。 CPUの中またはCPUパッケージの近くにはcacheがあります。 CPUは、memory controller を通してDRAMとデータをやり取りします。 一方、NVMe SSDは通常PCI Express (PCIe) に接続され、 SSD内部のNVMe controllerがNAND flashを制御します。

direction: right

cpu: {
  label: "CPU package"

  cores: {
    label: "CPU cores\ncore 0, core 1, ..."
  }

  cache: {
    label: "L1/L2/L3 cache"
  }

  memctl: {
    label: "memory controller"
  }

  pcie: {
    label: "PCIe root complex"
  }

  cores -> cache
}

dram: {
  label: "DRAM\nmain memory"
}

pch: {
  label: "chipset / PCH\noptional path"
}

ssd: {
  label: "NVMe SSD"

  nvme: {
    label: "NVMe controller"
  }

  nand: {
    label: "NAND flash"
  }

  nvme -> nand
}

cpu.memctl -> dram: "memory channels"
cpu.pcie -> ssd.nvme: "PCIe lanes\nCPU-attached"
cpu.pcie -> pch: "PCIe / chipset link"
pch -> ssd.nvme: "PCIe lanes\nchipset-attached"

ここでは、CPU package、DRAM、NVMe SSDの役割の違いだけを押さえます。 計算中に頻繁に読み書きする作業場所はDRAMで、SSDはファイルを保存する ストレージとして扱います。

実行中のデータ

プログラムを実行すると、入力ファイルや実行ファイルの内容はメモリに読み込まれます。 CPUはメモリ上のデータを読み、演算し、結果をメモリへ書き戻します。

storage  ->  memory  ->  CPU  ->  memory  ->  storage
 input       arrays      calc     results     output

保存されたCSV、HDF5、NumPy形式のファイルはストレージ上にあります。 計算中の Vec<f64> や ndarray::Array2<f64> はメモリ上にあります。 実際の演算はCPUで行われます。

計算量だけでは足りない

アルゴリズムを考えるとき、まず計算量を見ます。例えば、長さ n の配列の和は おおよそ O(n) の計算です。しかし、実行時間を考えるときは、 次のような要素も効きます。

• メモリから何バイト読むか。
• メモリへ何バイト書くか。
• データを連続して読めるか。
• 同じデータを何度も再利用できるか。
• ファイル入出力が計算時間に混ざっていないか。

このため、本書では計算本体、結果保存、plot を分けて考えます。 計算の速さを見たいときに、ファイル出力やplotの時間が混ざると、 何を測っているのかが分かりにくくなります。

latency・bandwidth・cache line

メモリアクセスには、主に2つの見方があります。

• latency: 1回の読み書きが始まってからデータが届くまでの待ち時間。
• bandwidth: 一定時間にどれだけ多くのデータを運べるか。

小さいデータを少しだけ読む場合は latency が効きやすく、 大きな配列を連続して読む場合は bandwidth が効きやすくなります。

cache

CPUはメモリよりはるかに高速に演算できます。 その差を埋めるために、CPUの近くには cache と呼ばれる小さく高速な記憶領域があります。 よく使うデータや、近くのデータをcacheに置くことで、 CPUがメモリを待つ時間を減らします。

重要なのは、CPUがメモリから1個の f64 だけを読むとは限らないことです。 多くの計算機では、メモリは cache line と呼ばれるまとまった単位でcacheへ運ばれます。 典型的なcache lineは64バイト程度です。

f64 は8バイトなので、64バイトのcache lineには f64 が8個入ります。 つまり、a[i] を読んだとき、近くの a[i + 1], a[i + 2] も 同じcache lineに入っている可能性があります。

連続アクセスが有利な理由

配列を先頭から順に読む場合、次に必要なデータが同じcache lineか、 次のcache lineにある可能性が高くなります。 このようなアクセスは計算機にとって予測しやすく、効率的です。

一方、遠く離れた要素を飛び飛びに読む場合、cache lineで運ばれたデータの多くを 使わないまま捨てることになります。

連続アクセス:
  a[0], a[1], a[2], a[3], ...

飛び飛びアクセス:
  a[0], a[1024], a[2048], ...

数値計算では、同じ数式でも配列の並びとloop orderによって速度が変わります。 これは、数学的な演算回数だけでなく、メモリからデータをどう運ぶかが効くためです。

compute bound と memory bandwidth bound

実行時間の主な制約が浮動小数点演算である場合、その計算は compute bound です。 一方、主な制約がメモリからデータを運ぶ速度である場合、 その計算は memory bandwidth bound です。

大きな配列に対して単純な演算を1回ずつ行う処理は、memory bandwidth bound になりがちです。 同じデータをcache上で何度も再利用できる処理は、compute bound に近づくことがあります。

連続アクセスと1次元データ

数値計算では、長い1次元配列を順に読む処理が頻繁に出てきます。 同じ要素数を読む場合でも、メモリ上で連続して読むか、飛び飛びに読むかで 実行時間が変わることがあります。

論理メモリ番地は1次元

プログラムから見るメモリ番地は、基本的に1次元の並びです。 プログラムは通常、連続した論理アドレス空間の中にデータが置かれているように扱えます。

address:  ... 1000 1008 1016 1024 1032 1040 ...
value:        a[0] a[1] a[2] a[3] a[4] a[5]

Vec<f64> の要素は、この1次元のアドレス空間の中で連続して並びます。 f64 は8バイトなので、隣の要素へ進むことは、典型的にはアドレスを8バイト進めることに対応します。

連続アクセス

次のように Vec<f64> を先頭から順に読むと、メモリ上でも隣の要素へ進みます。

fn sum_contiguous(x: &[f64]) -> f64 {
    let mut s = 0.0;
    for value in x {
        s += *value;
    }
    s
}

CPUはメモリから1個の f64 だけを読むのではなく、周辺のデータをまとめて cache line に載せます。 連続して読むと、cache line に載ったデータを無駄なく使いやすくなります。

飛び飛びアクセス

一方、次のように一定間隔で要素を読むと、実際に使わないデータも cache line に載ることがあります。

fn sum_every_k(x: &[f64], k: usize) -> f64 {
    let mut s = 0.0;
    let mut i = 0;
    while i < x.len() {
        s += x[i];
        i += k;
    }
    s
}

このようなアクセスでは、読んだバイト数のわりに使う値が少なくなり、 メモリ帯域やcacheの効率が悪くなることがあります。 ここでいう「何要素ごとに読むか」が、1次元データでのstrideの直感です。

実際には

実際には、論理アドレスはOSやMMUによってページ単位で物理メモリに対応づけられます。 ただし、配列をどの順序で読むかを考える段階では、まず論理アドレス空間での 連続性を見れば十分です。

第4章への接続

多次元配列では、数学的な添字が2個以上あっても、実データは1次元bufferに置かれます。 row-major、column-major、stride、loop order は、 多次元配列とメモリレイアウト で扱います。

stackとheap

実行中のプログラムが使うメモリには、大まかに stack と heap があります。 ここでは厳密な実装詳細ではなく、数値計算コードを読むための直感を押さえます。

stack

stack は、関数呼び出しに伴って使われる領域です。 小さい固定サイズの値や、関数の局所変数の一部が置かれます。 関数を抜けると、その関数のために使われていたstack領域はまとめて解放されます。 一方で、stack の大きさには上限があります。 Linux の一般的な環境では ulimit -s が 8192 KB、つまり約8 MiB に 設定されていることがよくあります。ただし、この値はOS、distribution、 shell の resource limit、thread の設定によって変わります。

例えば、次のような値は小さく、固定サイズです。

let dt: f64 = 0.01;
let n_steps: usize = 1000;
let x: [f64; 3] = [0.0, 1.0, 2.0];

数値計算では、巨大な固定長配列を局所変数として作ることは避けます。 例えば、次のような配列は f64 が100万個なので約8 MBあります。 これは stack の上限に近く、環境によっては stack overflow の原因になります。

let values: [f64; 1_000_000] = [0.0; 1_000_000];

大きな配列は、実行時に動的に確保し、heap に置くのが基本です。

heap

heap は、実行時にサイズが決まるデータや、大きなデータを置くための領域です。 数値計算で使う大きな配列は、通常heapに置かれます。

let n = 1_000_000;
let values = vec![0.0_f64; n];

Vec<f64> そのものは、長さ、容量、heap上のデータへのポインタを持つ小さな値です。 実際の100万個の f64 はheap側に置かれます。

所有権との関係

Rustでは、heap上のデータを誰が所有しているかを型システムで追跡します。 これは、後の章で扱う Vec<f64>、slice、借用の理解につながります。

• Vec<f64> はheap上のデータを所有する。
• &[f64] は所有せずに読む。
• &mut [f64] は所有せずに更新する。

この区別はRust固有の文法として現れますが、背景には 「大きなデータをどこに置き、誰が変更できるか」という一般的な問題があります。

第3章では、この考え方をRustの関数境界として扱います。

Rustで数値計算を書く最小セット

Rustと計算物理学では、 本書でRustを使う理由、開発環境、結果保存と可視化、本書の進め方を確認しました。 計算機の基本モデルでは、 CPU、メモリ、cache、stack、heap といった計算機の基本モデルを確認しました。 本章では、それをRustの数値計算コードとして書くための最小セットを整理します。

本章はRust Bookを置き換えるものではありません。対象を、計算物理のコードで繰り返し使う要素に絞ります。

• f64、usize、浮動小数点数
• 複素数と num-complex
• 関数、trait、let、mut、if、for
• Vec<f64>、配列、slice
• 所有権、借用、可変借用
• struct、Result
• 小さい単体テスト
• src/lib.rs、src/main.rs、module境界

AI agent を使う場合でも、これらの境界を理解しておくことは重要です。 agent に実装を任せる場合でも、関数の入力と出力、どのデータを所有するのか、 どこを可変にするのか、どの条件をテストするのかは、人間が確認する必要があります。

本章の構成

• 数値型・関数・小さいテスト f64、usize、関数、trait、制御構造、cargo testによる小さい検証を扱います。

• 浮動小数点演算と誤差 丸め誤差、桁落ち、情報落ち、浮動小数点数の比較、複素数、 Kahan summation などを扱います。

• Vecとスライス 1次元の数値データをVec<f64>で持ち、関数境界では&[f64]や&mut [f64]を使う方法を確認します。

• Resultとmodule境界 入力不正、配列サイズ不一致、収束失敗をResultで扱い、計算本体、入出力、実行用binaryを分ける方針を説明します。

本章を終えると、以降の章で出てくる小さな数値計算関数を読み、テストし、AI agent に安全に修正を任せるための基礎が整います。

数値型・関数・小さいテスト

数値計算では、まず小さな関数を正しく書き、それを手計算できる入力で検証することが重要です。本節では、以降の章で繰り返し使う最小限のRustの形を確認します。

なぜ小さく分けてテストするか

数値計算プログラムは、1回だけ数値を出せば正しい、というものではありません。 間違いを見つけ、原因を局所化し、修正できる形にしておく必要があります。

これは車を組み立てる作業に似ています。 車にはエンジン、ブレーキ、車輪、配線など多くの部品があります。 すべてを組み上げてから最後に走行テストだけを行うと、 不具合があることは分かっても、どの部品が原因かはすぐには分かりません。 実際には、部品ごとに検査し、いくつかの部品を組み合わせた段階でも検査し、 最後に全体を検査します。

ソフトウェアでも同じです。 小さい関数に分けて単体テストを書けば、計算の部品を個別に検査できます。 他の言語へ移植する場合も、プログラム全体を一度に移すのではなく、 部品ごとに移植し、同じ入力と期待値で正しさを確認できます。

テストを最初から導入することには、もう1つ重要な効果があります。 テストしやすい形にするためには、入力と出力が明確な関数に分ける必要があります。 つまり、早い段階でテストを書くと、自然に部品の境界が分かれます。 逆に、長いmain関数や隠れたグローバル状態に依存する形で書いた後から テストを足すのは困難です。

数値型

本書では、特に断りがない限り実数にはf64を使います。配列の添字や要素数にはusizeを使います。

let dt: f64 = 0.01;
let n_steps: usize = 1000;
let t = n_steps as f64 * dt;

浮動小数点数は、数学的な実数そのものではなく、有限のビット数で実数を近似した値です。 丸め誤差や比較の注意点については、 次節「浮動小数点演算と誤差」で扱います。

整数と浮動小数点数は自動では混ざりません。必要な場所で明示的に変換します。これは冗長に見えることもありますが、数値計算では、添字なのか物理量なのかを読み分けやすくする効果があります。

関数として切り出す

数値計算コードは、最初から大きなmain関数に書くのではなく、小さい関数に分けます。例えば、自由落下の位置を計算する関数は次のように書けます。

fn height(t: f64, h0: f64, v0: f64, g: f64) -> f64 {
    h0 + v0 * t - 0.5 * g * t * t
}

この関数は、入力と出力が明確です。ファイル入出力、プロット、コマンドライン引数とは独立しているため、単体テストしやすくなります。

traitの最小限の見方

trait（トレイト） は、ある型が「何をできるか」を表す約束です。 本書では、最初から複雑なジェネリクスを書く必要はありません。 まずは f64 や Vec<f64> のような具体的な型で関数を書き、必要になったときに trait を読みます。

数値計算でよく見る trait には、次のようなものがあります。

• Debug: {:?} で表示できる。
• Copy: 代入や関数呼び出しで値をコピーできる。
• Add, Sub, Mul, Div: +, -, *, / の演算に対応する。
• 外部crateのtrait: そのcrateが提供するメソッドを型に追加する。

例えば、後の章で出てくる ndarray-linalg では、 Norm や Solve のような trait を use することで、 配列にノルム計算や連立方程式を解くメソッドが見えるようになります。 このような use は、単に名前を短くするだけでなく、 trait が提供するメソッドを使うために必要な場合があります。

標準ライブラリには、f32 と f64 をまとめる Float trait はありません。 この教材では、特に理由がなければ f64 で書きます。 複数の浮動小数点型に対応する必要が出てきた場合だけ、num-traits のような 外部crateを検討します。

次の例は、trait bound を使った小さい関数です。 このような書き方は、関数を複数の数値型で使いたいときにだけ検討します。

use std::ops::{Add, Mul};

fn squared_norm2<T>(x: T, y: T) -> T
where
    T: Copy + Add<Output = T> + Mul<Output = T>,
{
    x * x + y * y
}

ここで T: Copy + Add<Output = T> + Mul<Output = T> は、 T がコピーでき、足し算と掛け算の結果も同じ T になることを要求しています。 ただし、最初からすべてを generic にする必要はありません。 物理量や配列の型が決まっている場合は、具体的な f64 の関数の方が読みやすいことも多いです。

小さいテスト

浮動小数点数では、==による厳密比較を避け、許容誤差を使って比較します。

fn height(t: f64, h0: f64, v0: f64, g: f64) -> f64 {
    h0 + v0 * t - 0.5 * g * t * t
}

fn close(a: f64, b: f64, tol: f64) -> bool {
    (a - b).abs() < tol
}

#[test]
fn height_at_initial_time_is_initial_height() {
    let h = height(0.0, 100.0, 0.0, 9.8);
    assert!(close(h, 100.0, 1e-12));
}

テストは、まず手計算できる小さい入力から始めます。解析解、保存量、極限ケース、対称性などがある場合は、それをテストに使います。

実行と検査

テストは次のコマンドで実行します。

cargo test

AI agent に実装を任せる場合も、少なくとも次の確認を行います。

cargo check
cargo test
git diff
git diff HEAD

cargo checkやcargo testは、型、借用、単体テストを検査します。 git diffはまだstageしていない変更を表示し、 git diff HEADは最終commitから現在の作業ツリーまでの変更全体を表示します。 ただし、コンパイルが通ることは、物理モデルや数値計算法が正しいことを意味しません。 人間は、問題設定、単位、境界条件、検証方法を確認します。

浮動小数点演算と誤差

計算物理学では、連続的な物理量をコンピュータ上で離散的に扱います。その際、実数を表現するために浮動小数点数が用いられますが、この表現は有限のビット数で行われるため、本質的に近似であり、誤差の発生は避けられません。

本節では、この浮動小数点数に起因する誤差について解説し、それが数値計算の結果に与える影響と、その対策について学びます。誤差の存在を正しく理解し、適切に対処することは、信頼性の高いシミュレーションを行う上で極めて重要です。

Rustにおける浮動小数点数

Rustの標準的な浮動小数点数型には、f32とf64の2種類があります。 どちらも数学的な実数そのものではなく、有限のビット数で実数を近似した値です。

f32はf64よりメモリ使用量が少ないため、大きな配列ではメモリ使用量や メモリ帯域の面で有利になることがあります。ただし、常に高速とは限りません。 速度は計算内容、ハードウェア、メモリアクセスのパターンに依存します。

Important

本書ではf64を標準とします

計算物理学の多くのシミュレーションでは、計算過程で発生する微小な誤差が積み重なり、最終的な結果に大きな影響を与えることがあります。そのため、本書では特に断りがない限り、より高い精度を持つf64型を標準として使用します。

複素数

Fourier 変換や量子力学では、複素数を使います。 Rustの標準ライブラリには複素数型が含まれていないため、 本書では外部クレート num-complex を使います。

Cargo.toml に以下を追加します。

[dependencies]
num-complex = "0.4"

num-complex では、実部と虚部が f64 の複素数型として Complex64 を使えます。

use num_complex::Complex64;

fn main() {
    let z = Complex64::new(1.0, 2.0); // 1 + 2i
    let w = Complex64::new(0.5, -1.0);

    println!("z + w = {}", z + w);
    println!("|z| = {}", z.norm());
}

複素数でも、実部と虚部は浮動小数点数です。 そのため、丸め誤差や比較の注意点は f64 と同じように現れます。 複素数の詳しい使い方は、Fourier 解析や量子力学の章で扱います。

浮動小数点数に起因する誤差

浮動小数点数を用いる際に発生する代表的な誤差を3つ紹介します。

1. 丸め誤差 (Rounding Error)

コンピュータが2進法で数値を扱うことに起因する誤差です。10進法では有限桁で表現できる小数（例: 0.1）が、2進法では無限小数となり、有限のビット数で正確に表現できない場合があります。その結果、最も近い表現可能な値に「丸め」られることで生じるのが丸め誤差です。

この影響を確認するために、簡単なコードを実行してみましょう。

fn main() {
    let a: f64 = 0.1;
    let b: f64 = 0.2;
    let sum = a + b;
    let expected: f64 = 0.3;

    println!("0.1 + 0.2 = {}", sum);
    println!("期待値     = {}", expected);
    println!("両者は等しいか？ => {}", sum == expected);

    // 内部的な表現を確認
    println!("実際の値 (内部表現): {:.20}", sum);
}

このコードを実行すると、以下のような結果が得られます。

0.1 + 0.2 = 0.30000000000000004
期待値     = 0.3
両者は等しいか？ => false
実際の値 (内部表現): 0.30000000000000004441

0.1 + 0.2 の計算結果が、期待する 0.3 とは微妙に異なる値になっていることがわかります。これは、0.1 と 0.2 が2進数で正確に表現できず、それぞれが近似値として格納されているために発生します。

Warning

浮動小数点数の直接比較は避ける

丸め誤差のため、== 演算子による浮動小数点数の直接比較は予期せぬ結果を招く可能性があります。計算結果が特定の値と等しくなることを期待するようなロジックは避けるべきです。比較の方法については後述します。

2. 桁落ち (Cancellation Error)

桁落ちとは、値がほぼ等しい2つの数値の差を計算した際に、有効桁数が大幅に失われてしまう現象です。

例えば、二次方程式 の解の公式を考えます。

ここで、 の場合、 となります。もし であれば、一方の解 の分子において、 と がほぼ等しい値となり、桁落ちが発生します。

以下のコードは、 の解をナイーブな解の公式で計算する例です。

fn main() {
    let a: f64 = 1.0;
    let b: f64 = 1.0e8; // 1億 (大きな値にすることで桁落ちが顕著になる)
    let c: f64 = 1.0;

    let d = (b * b - 4.0 * a * c).sqrt();

    // -b と d は非常に近い値になる (絶対値が近い)
    println!("-b = {}", -b);
    println!(" d = {}", d);

    // 桁落ちが発生する計算
    let x1_naive = (-b + d) / (2.0 * a);

    // もう一方の解
    let x2 = (-b - d) / (2.0 * a);

    println!("\nナイーブな計算:");
    println!("x1 = {:.17}", x1_naive);
    println!("x2 = {:.17}", x2);

    // 桁落ちを回避する工夫
    // 解と係数の関係 x1 * x2 = c/a を利用する
    let x1_stable = c / (a * x2);

    println!("\n安定な計算:");
    println!("x1 = {:.17}", x1_stable);

    // 相対誤差を計算
    let relative_error = (x1_naive - x1_stable).abs() / x1_stable.abs();
    println!(
        "\n相対誤差: {:.2e} ({:.1}%)",
        relative_error,
        relative_error * 100.0
    );
}

実行結果：

-b = -100000000
 d = 99999999.99999999

ナイーブな計算:
x1 = -0.00000000745058060
x2 = -100000000.00000000000000000

安定な計算:
x1 = -0.00000001000000000

相対誤差: 2.55e-1 (25.5%)

ナイーブな計算で得られた x1 は、本来の値（約 ）と比較して精度が大幅に低下しています。一方、解と係数の関係 を用いて と計算することで、桁落ちを回避し、より精度の高い解を得ることができました。

Note

なぜ x2 は桁落ちしないのか

x2 = (-b - d) / (2a) において、 のとき は負、 は正であるため、 は「負の値からさらに正の値を引く」計算となり、絶対値が増加する方向の演算です。桁落ちは「ほぼ等しい値の減算」で発生するため、x2 の計算では問題が生じません。このように、減算を回避する数式の変形が桁落ち対策の基本となります。

このように、計算順序やアルゴリズムを工夫することで、桁落ちを回避できる場合があります。

3. 情報落ち (Loss of Significance)

絶対値が非常に大きい数値と非常に小さい数値を足し算する際に、小さい数値が結果に反映されず、情報が失われてしまう現象を情報落ちと呼びます。

fn main() {
    let large_number: f64 = 1.0e16;
    let small_number: f64 = 1.0;

    let result = large_number + small_number - large_number;

    println!("(1.0e16 + 1.0) - 1.0e16 = {}", result);
    println!("期待値 = 1.0");
}

このプログラムの実行結果は 0.0 となります。1.0e16 + 1.0 の計算において、1.0 は f64 が表現できる有効桁数の範囲外となり、無視されてしまったためです。

多数の数値の和を計算する場合、絶対値の小さいものから順に足し合わせることで、情報落ちの影響を軽減できることがあります。より洗練された手法として、後述する Kahan summation（補正加算） があります。

特殊値：NaNとInfinity

浮動小数点数には、通常の数値以外に特殊な値が存在します。これらはシミュレーション中に意図せず発生することがあり、その挙動を理解しておくことが重要です。

Infinity（無限大）

f64::INFINITY（正の無限大）と f64::NEG_INFINITY（負の無限大）は、表現可能な範囲を超えた値を示します。

• オーバーフロー（非常に大きな数の演算）で発生
• 正の数をゼロで除算した場合に発生

NaN（Not a Number）

f64::NAN は、数学的に定義できない演算結果を示します。

• 0.0 / 0.0（ゼロ除算の不定形）
• (-1.0_f64).sqrt()（負の数の平方根）
• f64::INFINITY - f64::INFINITY（無限大同士の減算）

fn main() {
    // Infinityの発生
    println!("1.0 / 0.0 = {}", 1.0_f64 / 0.0);           // inf
    println!("1.0e308 * 10.0 = {}", 1.0e308_f64 * 10.0); // inf (オーバーフロー)

    // NaNの発生
    println!("0.0 / 0.0 = {}", 0.0_f64 / 0.0);           // NaN
    println!("(-1.0).sqrt() = {}", (-1.0_f64).sqrt());   // NaN

    // NaNの特殊な性質
    let nan = f64::NAN;
    println!("\nNaN == NaN は {}", nan == nan);  // false（重要！）
    println!("NaN.is_nan() = {}", nan.is_nan()); // true

    // 有限性のチェック
    let inf = f64::INFINITY;
    println!("\n1.0.is_finite() = {}", 1.0_f64.is_finite());   // true
    println!("inf.is_finite() = {}", inf.is_finite());         // false
    println!("nan.is_finite() = {}", nan.is_finite());         // false
}

Warning

NaNの伝播に注意

NaNを含む演算の結果はすべてNaNになります。シミュレーション中に一度でもNaNが発生すると、以降のすべての計算結果が「汚染」されてしまいます。

let nan = f64::NAN;
println!("{}", nan + 1.0);   // NaN
println!("{}", nan * 0.0);   // NaN
println!("{}", nan.sin());   // NaN

結果が期待と異なる場合や、シミュレーションが発散した場合は、is_nan() や is_finite() を用いて中間値をチェックすることが有効なデバッグ手法です。

let nan = f64::NAN;
println!("{}", nan.max(0.0)); // 0（NaNではない！）
println!("{}", nan.min(0.0)); // 0（NaNではない！）

Note

オーバーフローとアンダーフロー

f64 が表現できる最大値（約 ）を超えるとInfinityになります。一方、最小の正規化数（約 ）より小さい正の値は、ゼロまたは精度が低下した「非正規化数」になります（アンダーフロー）。

非常に大きな値や小さな値を扱う計算では、対数スケールでの計算や、適切なスケーリング（無次元化）を検討してください。

誤差への対策

これらの誤差が避けられないものである以上、それらを念頭に置いたプログラミングが求められます。

浮動小数点数の比較

前述の通り、== を用いた直接比較は危険です。代わりに、2つの数値の差がごく小さい値（許容誤差、イプシロンとも呼ばれる）未満であるかどうかで判定するのが一般的です。

絶対誤差による比較

最も単純な方法は、2つの数値の差の絶対値が許容誤差未満かを確認する方法です。

fn main() {
    let a: f64 = 0.1 + 0.2;
    let b: f64 = 0.3;

    let tolerance = 1e-10; // 許容誤差

    // 絶対誤差による比較
    if (a - b).abs() < tolerance {
        println!("絶対誤差の範囲で a と b は等しい");
    } else {
        println!("絶対誤差の範囲で a と b は等しくない");
    }
}

相対誤差による比較

絶対誤差による比較は、比較する数値が非常に大きい場合や非常に小さい場合にうまく機能しないことがあります。例えば、1.0e10 と 1.0e10 + 1.0 を比較する場合、その差 1.0 は 1e-10 よりもはるかに大きいですが、値のスケールから考えれば両者は「ほぼ等しい」と見なしたいかもしれません。

このような場合には、相対誤差を用いるのが有効です。

fn relative_eq(a: f64, b: f64, tolerance: f64) -> bool {
    // 注意: NaNが渡された場合、この関数はfalseを返す
    // （NaNとの比較はすべてfalseになるため）

    // ゼロ除算を避ける
    if a == 0.0 || b == 0.0 {
        return (a - b).abs() < tolerance;
    }
    (a - b).abs() / a.abs().max(b.abs()) < tolerance
}

この方法では、差を値の大きさで正規化するため、スケールに依存しない比較が可能になります。

Note

マシンイプシロン f64::EPSILON

Rustの標準ライブラリには、f64::EPSILON という定数が用意されています。これは「1.0と、1.0より大きい次の表現可能なf64の値との差」と定義されており、マシンイプシロンと呼ばれます。 ただし、この値は非常に小さく（約 ）、比較対象の数値のスケールによっては実用的でない場合があります。通常は、問題の性質に応じて適切な許容誤差 tolerance を独自に設定する必要があります。

実践的な比較: approx クレート

多くの場合、上記のような比較関数を自前で実装するよりも、approx のような実績のあるクレートを利用するのが最善です。このクレートは、絶対誤差、相対誤差、およびULP（Unit in the Last Place）に基づいた、頑健で使いやすい比較関数を提供します。

Cargo.toml に以下を追加して approx を導入できます。

[dependencies]
approx = "0.5"

そして、assert_relative_eq! マクロなどを使ってテストを書くことができます。

use approx::assert_relative_eq;

fn main() {
    let a: f64 = 0.1 + 0.2;
    let b: f64 = 0.3;

    // a と b が相対誤差の範囲で等しいことを表明する
    assert_relative_eq!(a, b, max_relative = 1e-10);

    println!("approxクレートを使って a と b がほぼ等しいことを確認しました。");
}

シミュレーションコード内の条件分岐などで比較を行いたい場合は、relative_eq!マクロやabs_diff_eq!マクロがブール値を返すため便利です。

本書では、テストコードなどで浮動小数点数の比較を行う際に approx クレートを積極的に利用します。

補正加算（Kahan Summation）

多数の数値を足し合わせる際、単純に加算を繰り返すと情報落ちによる誤差が蓄積します。Kahan summation（補正加算、カハンの加算アルゴリズム）は、誤差を追跡・補正しながら加算を行うことで、この問題を軽減するアルゴリズムです。

/// Kahan summationによる総和計算
fn kahan_sum(values: &[f64]) -> f64 {
    let mut sum = 0.0;
    let mut compensation = 0.0; // 累積誤差の補正項

    for &value in values {
        let y = value - compensation;  // 補正を適用した値
        let t = sum + y;               // 一時的な合計
        compensation = (t - sum) - y;  // 丸め誤差を記録
        sum = t;
    }
    sum
}

fn main() {
    // 小さな値を大量に足し合わせる例
    let n = 1_000_000;
    let small_value = 0.1;
    let values: Vec<f64> = vec![small_value; n];

    // 単純な加算
    let naive_sum: f64 = values.iter().sum();

    // Kahan summation
    let kahan = kahan_sum(&values);

    // 理論値
    let expected = small_value * n as f64;

    println!("理論値:           {:.17}", expected);
    println!("単純な加算:       {:.17}", naive_sum);
    println!("Kahan summation: {:.17}", kahan);
    println!("\n単純な加算の誤差: {:.2e}", (naive_sum - expected).abs());
    println!("Kahanの誤差:      {:.2e}", (kahan - expected).abs());
}

実行結果（環境により若干異なる場合があります）：

理論値:           100000.00000000000000000
単純な加算:       100000.00000133288267534
Kahan summation: 100000.00000000000000000

単純な加算の誤差: 1.33e-6
Kahanの誤差:      0.00e0

Kahan summation は、数値積分やモンテカルロ法など、多数の値を累積する計算で特に有効です。本書の後の章でも、精度が重要な場面でこのアルゴリズムを活用します。

数値的に安定なアルゴリズムの選択

同じ問題を解くアルゴリズムでも、計算手順によって誤差の蓄積度合いが異なる場合があります。桁落ちの例で見たように、数式を適切に変形したり、より数値的に安定したアルゴリズムを選択したりすることが重要です。

高精度計算の利用

f64 の精度でも不十分な場合は、より多くのビットを使って数値を表現する「多倍長浮動小数点数」ライブラリを利用するという選択肢もあります。

これについては、付録E「高精度計算」で詳しく解説します。

まとめ

本節では、計算物理学における基本でありながら、非常に重要な浮動小数点数の誤差について学びました。

• 丸め誤差・桁落ち・情報落ちといった誤差は、浮動小数点数演算において本質的に避けられない。
• NaNやInfinityといった特殊値の発生に注意し、必要に応じて is_nan() や is_finite() でチェックする。
• == による浮動小数点数の直接比較は危険であり、許容誤差を用いた比較が必要である。
• Kahan summation などの手法を用いることで、多数の値の加算における誤差を軽減できる。
• 計算順序の工夫や、数値的に安定なアルゴリズムの選択によって、誤差の影響を軽減できる。

シミュレーションを行う際は、常にこれらの誤差が結果に与える影響を意識することが、科学的に信頼できる結論を導くための第一歩となります。

Vecとスライス

計算物理では、時系列データ、粒子の状態、格子上の場など、多数の数値をまとめて扱います。1次元の数値データでは、まずVec<f64>とsliceを使えることが重要です。

所有するデータとしてのVec<f64>

Vec<f64>は、実行時に長さが決まる数値列を所有します。例えば、時刻ごとのエネルギーを保存する場合は次のように使います。

let mut energy: Vec<f64> = Vec::with_capacity(1000);

for step in 0..1000 {
    let e = step as f64 * 0.01;
    energy.push(e);
}

要素数の上限が分かっている場合は、Vec::with_capacityで容量を確保しておくと、途中の再割り当てを減らせます。

関数境界としての&[f64]

関数に数値列を渡すときは、多くの場合Vec<f64>そのものではなく&[f64]を受け取ります。

fn mean(xs: &[f64]) -> Option<f64> {
    if xs.is_empty() {
        return None;
    }

    let sum: f64 = xs.iter().sum();
    Some(sum / xs.len() as f64)
}

&[f64]はデータを所有しません。配列、Vec<f64>、配列の一部などを同じ関数に渡せます。

let xs = vec![1.0, 2.0, 3.0];
let avg = mean(&xs);
let first_two = mean(&xs[0..2]);

更新を伴う関数境界としての&mut [f64]

データをその場で更新する関数では、&mut [f64]を使います。

fn scale(xs: &mut [f64], factor: f64) {
    for x in xs {
        *x *= factor;
    }
}

let mut xs = vec![1.0, 2.0, 3.0];
scale(&mut xs, 0.5);

可変にする範囲を関数境界で明示することで、どのデータが変更されうるかを読み取りやすくなります。これはAI agent が生成した差分をレビューするときにも役立ちます。

stack、heap、所有権

小さな固定長配列、例えば3次元ベクトル[f64; 3]は値そのものとして扱いやすいデータです。一方、Vec<f64>は、長さが実行時に決まるデータをheap上に持ちます。

数値計算コードでは、次の使い分けを基本にします。

• 小さく固定長の量: [f64; 2], [f64; 3]
• 実行時に長さが決まる1次元データ: Vec<f64>
• 関数の入力: &[f64]
• 関数内で更新する入力: &mut [f64]

所有権や借用は、最初は複雑に見えるかもしれません。しかし、数値計算では「どの関数がデータを所有するか」「どの関数が読むだけか」「どの関数が更新するか」をはっきりさせる仕組みとして有用です。

Resultとmodule境界

数値計算コードでは、計算そのものだけでなく、入力不正、配列サイズ不一致、収束失敗、ファイル出力失敗なども扱う必要があります。小さい例ではpanic!でも動きますが、再利用する関数ではResultで失敗を返す方が安全です。

Resultで失敗を返す

例えば、2つのベクトルの内積では、長さが一致している必要があります。

fn dot(a: &[f64], b: &[f64]) -> Result<f64, String> {
    if a.len() != b.len() {
        return Err(format!(
            "length mismatch: a has {}, b has {}",
            a.len(),
            b.len()
        ));
    }

    let value = a.iter().zip(b).map(|(x, y)| x * y).sum();
    Ok(value)
}

呼び出し側では、失敗した場合に何をするかを明示できます。

let a = vec![1.0, 2.0, 3.0];
let b = vec![4.0, 5.0, 6.0];
let value = dot(&a, &b)?;

structでパラメータをまとめる

引数が増えてきたら、パラメータや状態をstructにまとめます。

struct EulerConfig {
    dt: f64,
    n_steps: usize,
}

impl EulerConfig {
    fn validate(&self) -> Result<(), String> {
        if self.dt <= 0.0 {
            return Err("dt must be positive".to_string());
        }
        if self.n_steps == 0 {
            return Err("n_steps must be positive".to_string());
        }
        Ok(())
    }
}

入力条件をvalidateに分けると、テストしやすくなります。

module境界

大きい計算では、main.rsにすべてを書かないようにします。例えば、次のように分けます。

src/
├── lib.rs
├── model.rs
├── algorithm.rs
├── io.rs
└── main.rs

役割は次のように分けます。

• model: 物理モデル、パラメータ、状態量
• algorithm: 時間発展、反復法、更新規則
• io: 結果保存、metadata保存
• main: コマンドライン引数、実行順序
• tests: 小さい検証問題、解析解との比較

これは、前節で述べた「小さい部品を先に検査する」考え方を、 project全体に広げたものです。 model、algorithm、ioを分けておくと、各部品を単体テストしやすくなります。 他の言語へ移植する場合も、moduleごとに移し、同じ小さい検証問題で比較できます。

この分割は、AI agent に実装を任せる場合にも有効です。「algorithmだけを修正してください」「ioにCSV保存を追加してください」「modelの物理量更新を疑似コード化してください」のように、作業範囲を限定できます。

panic!してよい場合

panic!は、プログラムのバグを示す場合や、サンプルコードを短くする場合には使えます。一方、ユーザー入力、ファイル入出力、収束失敗、配列サイズ不一致のように実行時に起こりうる失敗は、Resultで返す方が適しています。

本書では、説明を短くするためにexpectやpanic!を使うこともあります。ただし、研究用・演習用のコードを長く育てる場合は、どの失敗を呼び出し側に返すべきかを意識してください。

多次元データと配列

Important

本章の前提知識

本章の内容を理解するためには、以下の章を事前に学習しておく必要があります。

• Rustと計算物理学
• 計算機の基本モデル
• Rustで数値計算を書く最小セット

本章では、1次元の数値データから、多次元配列、memory layout、外部クレート、 データ保存へ進みます。以降の線形代数、偏微分方程式、統計力学、量子力学では、 多次元データを正しく扱うことが重要になります。

本章の目的は、ndarrayを使うことだけではありません。 計算機の基本モデルで見た cache line、連続アクセス、stride の考え方を使い、 配列のshape、axis、copy/view、保存形式、metadataを意識して、 計算結果を後から検証できる形で扱うことを目指します。

本章の構成

• 1次元データから多次元データへ Rust標準の配列、スライス、Vec<T>を出発点として、ベクトル、行列、格子データへ進む準備をします。

• 多次元配列とメモリレイアウト flattening、row-major、column-major、stride、view/copyの違いをRustの配列で確認します。

• 外部クレートの活用（ndarray入門） Rustの科学技術計算で標準的に使われる軽量なndarrayクレートと、 LAPACK系の線形代数を使うndarray-linalgの位置づけを確認します。

• 発展的なテンソルライブラリ tenferro ndarrayを確認した後、自動微分、GPU実行、テンソル縮約を扱う 発展的な選択肢としてtenferroを確認します。

• データ保存とmetadata CSV、HDF5、NumPy形式、metadataの使い分けを整理します。

本章を終えることで、読者は多次元の数値データをRustで扱い、保存し、後続の解析や可視化につなげるための基礎を習得できます。

1次元データから多次元データへ

Note

本節のポイント

• Rust標準の配列 ([T; N])、スライス (&[T])、ベクタ (Vec<T>) を多次元データの入口として整理する。
• 1次元の連続データから、ベクトル、行列、格子データへ進む考え方を確認する。
• 次節で扱うmemory layoutやndarrayの必要性を理解する。

Rustのこれらコレクション型に習熟している方は、本節を読み飛ばして次節の「多次元配列とメモリレイアウト」に進んでも構いません。

物理シミュレーションやデータ解析では、個々の数値を扱うだけでなく、多数の数値をひとまとめにしてベクトルや行列として扱う場面が頻繁に現れます。例えば、3次元空間における粒子の位置や速度、あるいは観測データの時系列などがそれに当たります。

Rustの標準ライブラリは、このようなデータの集合を扱うための基本的なツールを提供しています。本節では、その中でも特に重要な配列 ([T; N])、スライス (&[T])、そしてベクタ (Vec<T>) の3つを取り上げ、その特徴と基本的な使い方、そして数値計算における使い分けについて解説します。

これらの基本的なデータ構造を理解することは、より高度な科学技術計算ライブラリ（例えば後の節で学ぶndarray）を効果的に利用するための土台となります。

配列 [T; N]

配列は、最も基本的なコレクション型です。[T; N]という型で表され、Tは要素の型、Nは要素数を意味します。

特徴

• 固定長: 配列の長さ（要素数）はコンパイル時に決定され、後から変更することはできません。
• スタック確保: 配列のデータは、通常スタック上にメモリが確保されます。これにより、ヒープ割り当てのオーバーヘッドがなく、非常に高速にアクセスできます。

数値計算の文脈では、3次元ベクトルや4元数（クォータニオン）のように、コンパイル時に次元数が確定している小さなデータ構造を表現するのに適しています。

基本的な使い方

配列の宣言と初期化は、以下のように行います。

// 3次元の位置ベクトルを表す配列
let position: [f64; 3] = [1.0, 2.0, 3.0];

// すべての要素を0.0で初期化
let velocity: [f64; 3] = [0.0; 3];

// 要素へのアクセス
println!("x-coordinate: {}", position[0]);
println!("Initial velocity: {:?}", velocity);

// forループによるイテレーション
for component in &position {
    println!("{}", component);
}

長所と短所

• 長所: 高速なアクセス。コンパイル時のサイズチェックによる安全性。
• 短所: サイズが固定されているため、実行時に要素数を変更したい場合には不向きです。

スライス &[T]

スライスは、それ自体がデータを所有するのではなく、配列やVecといった他のコレクション型の一部または全体を「参照」するための型です。この「参照する」という性質がスライスの最も重要な特徴です。

特徴

• 所有権を持たない (Borrowing): スライスはデータの所有権を持たず、データを借用（borrow）します。これにより、大きなデータブロックをコピーすることなく、関数にデータの部分集合を効率的に渡すことができます。
• 動的な長さ: スライスが参照する範囲の長さはコンパイル時には決まっておらず、実行時に決まります。

スライスは、あるデータ集合に対して操作を行う関数を定義する際に、非常に柔軟で効率的なインターフェースを提供します。

基本的な使い方

スライスは、配列やVecから簡単に作成できます。

fn process_data(data: &[f64]) {
    // dataスライス内のすべての要素の平均値を計算する
    let sum: f64 = data.iter().sum();
    let count = data.len();
    if count > 0 {
        println!("Average: {}", sum / count as f64);
    }
}

fn main() {
    let array: [f64; 5] = [1.0, 2.0, 3.0, 4.0, 5.0];
    let vector: Vec<f64> = vec![10.0, 20.0, 30.0];

    // 配列全体をスライスとして渡す
    process_data(&array);

    // 配列の一部をスライスとして渡す
    process_data(&array[1..4]); // インデックス1から3までの要素 ([2.0, 3.0, 4.0])

    // Vec全体をスライスとして渡す
    process_data(&vector);
}

process_data関数は、引数として&[f64]（f64のスライス）を取ります。これにより、この関数は配列からもVecからも、その一部からでもデータを受け取ることができます。これは非常に強力な抽象化です。

また、&mut [T]を用いることで、スライスを通じて元のデータを変更することも可能です。

fn scale(data: &mut [f64], factor: f64) {
    for element in data {
        *element *= factor;
    }
}

fn main() {
    let mut values = vec![1.0, 2.0, 3.0];
    scale(&mut values, 2.0);
    println!("{:?}", values); // [2.0, 4.0, 6.0]
}

長所と短所

• 長所: 配列とVecの両方を同じインターフェースで扱える汎用性。データをコピーせずに部分参照が可能。
• 短所: 参照先のライフタイムを考慮する必要がある。所有権を持たないため、データを長期間保持する用途には使えない。

ベクタ Vec<T>

ベクタは、可変長のコレクション型です。データの実体はヒープ上に確保され、実行時に要素を追加したり削除したりすることができます。

特徴

• 可変長: pushメソッドで要素を追加したり、popメソッドで末尾の要素を削除したりと、サイズを自由に変更できます。
• ヒープ確保: データがヒープ上に存在するため、コンパイル時にサイズが不明なデータを扱うことができます。

数値計算では、シミュレーションの各ステップで得られる結果を時系列データとして記録していく場合など、実行時に要素数が動的に変化する状況で非常に役立ちます。

基本的な使い方

// 空のベクタを作成
let mut time_series: Vec<f64> = Vec::new();

// vec!マクロで初期値を指定して作成
let positions = vec![[0.0, 0.0], [1.0, 1.0]];
println!("Initial positions: {:?}", positions);

// あらかじめ容量を確保しておくと再割り当てを回避できる
let results: Vec<f64> = Vec::with_capacity(1000);

// push: 要素を追加
time_series.push(0.0);
time_series.push(0.1);
time_series.push(0.2);

// 要素へのアクセス
println!("Latest data point: {}", time_series[2]);

// forループによるイテレーション
for data_point in &time_series {
    println!("{}", data_point);
}

println!("Number of data points: {}", time_series.len());
println!("Results capacity: {}", results.capacity());

Tip

境界チェックについて

Rustでは、配列やVecへのインデックスアクセス（例: time_series[0]）は実行時に境界チェックが行われます。範囲外のインデックスにアクセスしようとすると、プログラムはパニックして停止します。これは、C/C++で起こりうるバッファオーバーフローのような未定義動作を防ぐための安全機構です。

長所と短所

• 長所: 実行時にサイズを柔軟に変更できる高い柔軟性。
• 短所: ヒープ割り当てに伴うわずかなパフォーマンス上のオーバーヘッド。要素追加時に容量が不足すると、より大きなメモリ領域を確保してデータをコピーする「再割り当て（reallocation）」が発生する可能性があります。あらかじめ要素数の上限が分かっている場合は、Vec::with_capacityで容量を確保しておくことで再割り当てを回避できます。

まとめと使い分け

これら3つのデータ構造の特性と、数値計算における主な用途をまとめます。

これらの使い分けを理解することは、効率的で読みやすいコードを書くための第一歩です。

• 配列 [T; N]: 扱うデータの次元や要素数が初めから分かっている場合に最適です。
• Vec<T>: データの数が実行中に変わる可能性がある場合に選択します。
• スライス &[T]: 関数にデータの集合を渡す際の引数の型として、ほとんどの場合に最良の選択です。これにより、関数がより汎用的になり、配列とVecの両方で再利用できます。

次のステップへ

本節では、Rust標準ライブラリが提供する基本的なデータ構造について学びました。これらを用いることで、基本的なベクトル演算やデータ管理は可能です。

しかし、より複雑な科学技術計算、例えば行列、格子、場のデータを扱うには、これらのプリミティブな型だけでは表現が冗長になりやすく、shapeやaxisの管理も難しくなります。

そこで次節では、多次元の添字がメモリ上でどのように並ぶかを確認します。その上で、Rustの科学技術計算エコシステムにおける中核的なクレートであるndarrayへ進みます。

多次元配列とメモリレイアウト

格子、行列、画像、波動関数、相関関数などを扱うとき、データは多次元の添字を持ちます。 しかし、計算機の基本モデルで見たように、 コンピュータのメモリ上ではデータは基本的に1次元に並びます。 この対応を理解しておくと、配列の扱い、性能、ファイル保存を考えやすくなります。

flattening

2次元配列u[i, j]を1次元のVec<f64>に保存する場合、例えば行優先（row-major）では次のように対応させます。

fn idx(i: usize, j: usize, nx: usize) -> usize {
    i * nx + j
}

このとき、iは行、jは列、nxは1行あたりの要素数です。

let ny = 4;
let nx = 5;
let mut u = vec![0.0; ny * nx];

let i = 2;
let j = 3;
u[idx(i, j, nx)] = 1.0;

flattening を明示すると、メモリ上の並び、境界条件、ファイル出力の対応を確認しやすくなります。

row-major と column-major

多次元配列の並べ方には、主に次の2つがあります。

• row-major: 同じ行の要素が連続する。C/C++、ndarray、NumPyの標準に近い。
• column-major: 同じ列の要素が連続する。Fortran、MATLAB、Julia、LAPACK/BLAS、Eigen3でよく使われる。

どちらが正しいという話ではありません。重要なのは、配列を作る側、計算する側、保存する側、読む側で同じ約束を使うことです。

ここでは、ny 行 nx 列の配列を考えます。 i を行index、j を列indexとすると、row-major では典型的に次の対応になります。

index(i, j) = i * nx + j

このとき、j を1つ増やすとメモリ上でも隣の要素へ進みます。 i を1つ増やすと nx 要素分だけ進みます。

column-major では、同じ形の配列を次のように対応させます。

index(i, j) = j * ny + i

このとき、i を1つ増やすとメモリ上でも隣の要素へ進みます。 j を1つ増やすと ny 要素分だけ進みます。

stride

stride は、ある軸に沿って添字を1つ進めたとき、メモリ上で何要素進むかを表します。 連続した配列ではアクセスが速く、stride が大きいアクセスでは cache line に載ったデータを十分に使えないことがあります。

ny 行 nx 列の2次元配列では、典型的には次のようになります。

row-major:

• j 方向の stride: 1
• i 方向の stride: nx

column-major:

• i 方向の stride: 1
• j 方向の stride: ny

例えば row-major の2次元配列では、同じ行の隣接要素を読むループは連続アクセスになります。

for i in 0..ny {
    for j in 0..nx {
        let value = u[idx(i, j, nx)];
        // valueを使う
    }
}

計算量だけでなく、メモリアクセスの順序も性能に影響します。 大きな配列を扱う章では、loop order、cache、memory bandwidth を意識します。

逆に、row-major の配列を列方向に読むと、nx 要素ずつ飛ぶアクセスになります。

for j in 0..nx {
    for i in 0..ny {
        let value = u[idx(i, j, nx)];
        // valueを使う
    }
}

小さい配列では差が見えないこともあります。 しかし、大きな配列や何度も繰り返す計算では、この違いが実行時間に効きます。

N次元の場合

多次元配列でも考え方は同じです。 shape を (n0, n1, ..., n_{d-1})、添字を (i0, i1, ..., i_{d-1}) と書くと、index は各添字とstrideの積の和になります。

index = i0 * stride0 + i1 * stride1 + ... + i_{d-1} * stride_{d-1}

row-major では右端の添字が最も速く変わります。

strides = (n1 * n2 * ... * n_{d-1}, ..., n_{d-1}, 1)

column-major では左端の添字が最も速く変わります。

strides = (1, n0, n0 * n1, ..., n0 * n1 * ... * n_{d-2})

この規約は、数学的な shape そのものとは別です。 同じ 2 x 3 の行列でも、row-major と column-major では 1次元bufferへの並び方が変わります。 ndarray のような配列型では、shape に加えて stride などのmetadataを持つことで、 同じbufferをさまざまな見方で扱えます。

view、copy、reshape

多次元配列ライブラリでは、部分配列や転置を「view」として扱える場合があります。view は元データへの参照であり、データをコピーしません。一方、連続した新しい配列が必要な場合は、実データのコピーが発生します。

AI agent に配列操作を任せる場合は、次の点を確認してください。

• view で十分な場所で不要な copy をしていないか。
• reshape 後の shape とデータ順序が期待通りか。
• transpose 後に、後続の処理が想定する memory layout と合っているか。
• 保存時に shape、axis、単位が metadata として残っているか。

多次元配列のバグは、コンパイルは通っても、軸の取り違えや境界条件の間違いとして現れます。小さい配列で手計算できるテストを作ることが重要です。

外部クレートの活用（ndarray入門）

Note

本節のポイント

• Rustの科学技術計算における標準的で軽量な配列クレートndarrayの基本を学ぶ。
• ndarrayを用いて、ベクトルや行列といった多次元配列を効率的に作成・操作する方法を習得する。
• ndarray-linalgが、LAPACK系の線形代数をndarrayに接続するcrateであることを確認する。

本節は、前節「1次元データから多次元データへ」と「多次元配列とメモリレイアウト」で学んだ内容を土台としています。

前の節までで、Rustの標準ライブラリが提供する配列、スライス、Vec<T>、そして多次元データのmemory layoutを確認しました。これらはデータ管理の基本として有用ですが、より高度な科学技術計算、特に多次元配列（ベクトル、行列、テンソル）を多用する線形代数演算には力不足です。

例えば、Vec<Vec<f64>>で「行列」を表現しようとすると、以下のような課題に直面します。

• メモリ効率: 各行が独立したヒープ割り当てとなり、データがメモリ上で連続的に配置されないため、キャッシュ効率が悪化します。
• 演算の煩雑さ: 行列の積や転置といった基本的な演算も、手動でループを実装する必要があり、コードが冗長かつエラーが発生しやすくなります。
• パフォーマンス: 最適化された線形代数ライブラリ（BLASなど）の恩恵を受けにくく、手実装のループでは高いパフォーマンスを得るのが困難です。

このような課題を解決する基本的な選択肢が、ndarray クレートです。 ndarrayは、RustでN次元配列、array view、多次元slicing、効率的な演算を扱うための 標準的で軽量なcrateです。

さらに、固有値分解、特異値分解、線形方程式の求解などの既存LAPACK系の線形代数を ndarray と接続したい場合は、ndarray-linalg が候補になります。 これは ndarray の ArrayBase に線形代数機能を提供し、 外部のLAPACK実装を利用するcrateです。

本節では、まずndarrayの基本的な使い方を学び、 必要に応じてndarray-linalgへ進む位置づけを確認します。

ndarrayのセットアップ

まず、ndarrayをプロジェクトに追加します。Cargo.tomlファイルの[dependencies]セクションに以下の行を追記してください。

[dependencies]
ndarray = "0.17"

バージョンは執筆時点のものです。最新版はcrates.ioで確認できます。

ndarray-linalgを使う場合は、LAPACK backendも選ぶ必要があります。 OpenBLAS、Netlib、Intel MKLなどのbackend featureから1つを選ぶ形です。 教材では詳細なbackend構築手順には踏み込みません。 必要になった時点で、公式documentationを確認して、使用環境に合うbackendを選びます。

本書のコード例では、個別の型、関数、traitをuse宣言で明示的にインポートします。 これにより、どの機能がどのcrateから来ているかを追いやすくします。

use ndarray::{arr1, arr2, Array, Array1, Array2, Axis};

1次元配列 (Array1)：ベクトルの表現

ndarrayにおける1次元配列Array1<T>は、数学的なベクトルに対応します。

ベクトルの作成

ndarrayには、ベクトルを生成するための多様な方法が用意されています。

use ndarray::{arr1, Array, Array1};

fn main() {
    // スライスからベクトルを作成
    let a: Array1<f64> = Array::from(vec![1.0, 2.0, 3.0]);
    println!("a = {}", a);

    // arr1 関数を使った、より簡潔な作成方法
    let b = arr1(&[4.0, 5.0, 6.0]);
    println!("b = {}", b);

    // 0から9までの連番を持つベクトルを作成
    let c: Array1<f64> = Array::range(0.0, 10.0, 1.0);
    println!("c = {}", c);

    // すべての要素が0.0のベクトルを作成（次元を指定）
    let zeros = Array1::<f64>::zeros(5);
    println!("zeros = {}", zeros);
}

ベクトル演算

ndarrayの最大の利点の一つは、ベクトルや行列に対する数学的な演算が直感的に書ける点です。+, -, *, / といった演算子がオーバーロードされており、要素ごとの演算が可能です。

use ndarray::arr1;

fn main() {
    let a = arr1(&[1.0, 2.0, 3.0]);
    let b = arr1(&[4.0, 5.0, 6.0]);

    // ベクトルの加算
    println!("a + b = {}", &a + &b);

    // ベクトルの減算
    println!("a - b = {}", &a - &b);

    // スカラー倍
    println!("a * 2.0 = {}", &a * 2.0);

    // 要素ごとの積
    println!("a * b (element-wise) = {}", &a * &b);

    // 内積 (dot product)
    let dot_product = a.dot(&b);
    println!("a . b = {}", dot_product);
}

Note

参照と所有権 &a + &bのように、演算時にベクトルの参照(&)を渡している点に注目してください。ndarrayの多くの二項演算は参照を受け取るように実装されています。これにより、演算後も元のデータaとbの所有権が維持され、続けて使用することができます。もしa + bと記述すると、aの所有権が移動してしまい、以降aは使えなくなります。

2次元配列 (Array2)：行列の表現

ndarrayにおける2次元配列Array2<T>は、行列に対応します。

行列の作成

行列もベクトルと同様、様々な方法で作成できます。 ndarray のowned arrayは、標準では row-major、つまり同じ行の隣接要素が メモリ上で連続する配置です。これはC-orderとも呼ばれます。 Array2::from_shape_vec((rows, cols), data) に渡す data も、この順序で解釈されます。

use ndarray::{arr2, Array2};

fn main() {
    // arr2 関数で行列を作成
    // 各内部配列が「行」に対応する
    let m = arr2(&[[1.0, 2.0, 3.0],
                   [4.0, 5.0, 6.0]]);
    println!("m =\n{}", m);

    // from_shape_vecで行列を作成
    // (行数, 列数)のタプルと、データを格納したVecを渡す
    let shape = (3, 2);
    let data = vec![1, 2, 3, 4, 5, 6];
    let m_from_vec = Array2::from_shape_vec(shape, data)
        .expect("Incompatible shape");
    println!("m_from_vec =\n{}", m_from_vec);

    // 3x3のゼロ行列
    let zeros = Array2::<f64>::zeros((3, 3));
    println!("zeros =\n{}", zeros);

    // 2x2の単位行列
    let eye = Array2::<f64>::eye(2);
    println!("eye =\n{}", eye);
}

Note

layoutとstride

ndarray では、配列がshapeとstrideを持ちます。 owned arrayを普通に作る場合は row-major と考えてよいですが、 slice、transpose、view ではstrideが変わります。 性能や外部ライブラリとの接続が重要な場合は、shapeだけでなくstrideも確認します。

行列演算

行列同士の加減算やスカラー倍もベクトルと同様に直感的に記述できます。そして、dotメソッドを用いることで、行列積や行列とベクトルの積を計算できます。

use ndarray::{arr1, arr2};

fn main() {
    let m1 = arr2(&[[1.0, 2.0],
                    [3.0, 4.0]]);

    let m2 = arr2(&[[5.0, 6.0],
                    [7.0, 8.0]]);

    // 行列の加算
    println!("m1 + m2 =\n{}", &m1 + &m2);

    // 行列のスカラー倍
    println!("m1 * 3.0 =\n{}", &m1 * 3.0);

    // 行列積
    let matrix_product = m1.dot(&m2);
    println!("m1 * m2 (matrix product) =\n{}", matrix_product);

    // 行列とベクトルの積
    let v = arr1(&[10.0, 20.0]);
    let matrix_vector_product = m1.dot(&v);
    println!("m1 * v =\n{}", matrix_vector_product);
}

ndarray 自体でも、基本的な行列積や要素ごとの演算を扱えます。 一方、分解、固有値、特異値分解、線形方程式の求解など、 LAPACK系の機能が必要な場合は ndarray-linalg を使います。 ndarray-linalg はbackend featureと外部libraryの設定に依存するため、 小さい演習で不用意に必須化しません。

スライシング

ndarrayの特に強力な機能がスライシングです。これにより、配列の一部を効率的に、かつ柔軟に抜き出すことができます。スライシングはNumPyのそれに非常によく似ており、s![]マクロを用いて行います。

use ndarray::{arr2, s};

fn main() {
    let m = arr2(&[[ 1,  2,  3,  4],
                   [ 5,  6,  7,  8],
                   [ 9, 10, 11, 12]]);
    println!("m =\n{}", m);

    // 0番目の行を取得
    let row0 = m.slice(s![0, ..]);
    println!("Row 0: {}", row0); // [1, 2, 3, 4]

    // 1番目の列を取得
    let col1 = m.slice(s![.., 1]);
    println!("Column 1: {}", col1); // [2, 6, 10]

    // 部分行列を取得
    // 0-1行目、1-2列目を抜き出す
    let sub_matrix = m.slice(s![0..2, 1..3]);
    println!("Sub-matrix (0..2, 1..3) =\n{}", sub_matrix);
    // [[2, 3],
    //  [6, 7]]

    // スライスは元のデータのビュー（参照）
    // スライスを介して元のデータを変更することも可能
    let mut m_mut = m.clone();
    let mut sub_view = m_mut.slice_mut(s![0..2, 1..3]);
    sub_view.fill(0); // 部分行列を0で埋める
    println!("Modified m_mut =\n{}", m_mut);
}

Important

スライスはビューである

slice()メソッドが返すArrayView型は、元の配列データの所有権を持たない「ビュー（参照）」です。これにより、巨大な行列の一部を操作する際に、不要なデータコピーを回避し、高いメモリ効率とパフォーマンスを実現します。

ブロードキャスト

ブロードキャストは、形状（shape）が異なる配列間の演算を可能にする強力な機能です。ndarrayは、不足している次元を自動的に「引き伸ばし」て、形状を一致させてから演算を実行します。

例えば、行列のすべての要素に同じ値を加算する、あるいは行列の各行に同じベクトルを加算するといった操作が簡単に行えます。

use ndarray::{arr1, arr2, Axis};

fn main() {
    let m = arr2(&[[1.0, 2.0],
                   [3.0, 4.0]]);

    let v = arr1(&[10.0, 20.0]);

    // 行列の各行にベクトルvを加算する
    // vは (2,) -> (1, 2) にブロードキャストされ、
    // それがさらに (2, 2) に引き伸ばされてmと加算される
    let result_row = &m + &v;
    println!("m + v (row-wise broadcast) =\n{}", result_row);

    // 行列の各列にベクトルを加算する場合
    // ベクトルを列ベクトルとして扱うために次元を追加する必要がある
    // v.view().insert_axis(Axis(1)) は vの形状を (2,) から (2, 1) に変える
    // view() を使うことで v の所有権を保持したままビューを操作できる
    let result_col = &m + &v.view().insert_axis(Axis(1));
    println!("m + v (column-wise broadcast) =\n{}", result_col);
}

実行結果:

m + v (row-wise broadcast) =
[[11, 22],
 [13, 24]]
m + v (column-wise broadcast) =
[[11, 12],
 [23, 24]]

ndarrayとRust標準型の連携

ndarrayは標準ライブラリの型とスムーズに連携できます。

• ArrayからVecへ: to_vec()
• Arrayからスライスへ: as_slice()
• スライスからArrayViewへ: ArrayView::from()

これにより、既存のRustコードにndarrayを段階的に導入したり、ndarrayで計算した結果を他のライブラリに渡したりすることが容易になります。

その他の高度な機能

ndarrayは本節で紹介した以外にも、以下のような多数の機能を提供します。

• 高階関数: map, fold, zipなど、配列を効率的に操作するメソッド。
• 軸に沿った操作: sum, mean, max, minなどを特定の軸（行方向や列方向）に沿って計算。
• LAPACK系の線形代数: ndarray-linalgクレートと組み合わせることで、分解、固有値問題、特異値分解（SVD）、線形方程式の求解などを扱えます。

まとめ

本節では、Rustの科学技術計算でよく使われる軽量な配列クレートndarrayの基本的な使い方を学びました。

• Array1（ベクトル）とArray2（行列）の作成方法と、それらの間の基本的な演算（四則演算、内積、行列積）を習得しました。
• s![]マクロを用いた柔軟なスライシング機能により、配列の一部を効率的に参照・操作できることを見ました。
• ブロードキャスト機能により、形状の異なる配列間でも直感的な演算が可能であることを学びました。

ndarrayは、Rustで多次元配列を扱う標準的な軽量crateです。 LAPACK系の線形代数が必要な場合は ndarray-linalg も確認します。

参照:

• ndarray: https://github.com/rust-ndarray/ndarray
• ndarray-linalg: https://github.com/rust-ndarray/ndarray-linalg
• ndarray-linalg API: https://rust-ndarray.github.io/ndarray-linalg/ndarray_linalg/

発展的なテンソルライブラリ tenferro

Note

本節のポイント

• まずは Vec<f64>、&[f64]、ndarray を基本にする。
• tenferro は、より高機能なテンソル計算が必要になったときの発展的な選択肢である。
• tenferro は本書の著者の1人が開発に関わるライブラリである。
• tenferro は単一のクレートではなく、tenferro-tensor や tenferro-ad などの 複数クレートからなる workspace で、crates.io に公開されている。
• tenferro の密なテンソルは column-major（列優先）なので、ndarray の標準的な row-major（行優先）と混同しない。
• tenferro はまだ 0.x 版（pre-1.0）なので、使う場合は公式リポジトリとドキュメントで 現在のAPIを確認する。

1次元の所有データには、通常 Vec<f64> を使います。 関数に渡すときは、&[f64] や &mut [f64] を使うのが基本です。

2次元以上の配列を軽く扱うだけなら、ndarray が自然な選択肢です。 本書の多くの例でも、まず ndarray を使います。

tenferro は、それより発展的なテンソル計算のためのライブラリです。 ここでいうテンソルは、ベクトルや行列をさらに多次元にした数値配列だと思えば十分です。 単に2次元配列を持つためだけなら、tenferro は大きすぎる選択になりがちです。

本書で tenferro を紹介する理由の1つは、本書の著者の1人が開発に関わっているためです。 Rustで科学技術計算を発展させるときに、どのような機能が必要になり、 どこに実装上の注意があるかを考える具体例として扱います。

どのような場合に使うか

tenferro は、次のような機能を同じ枠組みで扱いたい場合に候補になります。

• 自動微分: 計算式から、微分や勾配を自動的に計算する機能です。 最適化問題や逆問題で使うことがあります。
• GPU実行: CPUではなくGPU上で大きなテンソル演算を実行します。 ただし、CPUとGPUの間でデータを移す時間も考える必要があります。
• テンソル縮約: 添字を持つ多次元配列について、ある添字に関して和を取る演算です。 行列積もテンソル縮約の一例です。 Python/NumPy では einsum という名前でよく使われます。
• 線形代数やFFTとの接続: 行列分解、連立方程式、フーリエ変換などをテンソル計算の流れの中で扱います。

このような機能が必要でなければ、まず ndarray と ndarray-linalg で十分かを考えます。

クレート構成

tenferro は単一の tenferro クレートではありません。 tenferro-tensor、tenferro-ad、tenferro-einsum、tenferro-linalg、 tenferro-fft、tenferro-cpu、tenferro-runtime などの複数クレートからなる workspace で、それぞれが crates.io に公開されています。 ルートの facade クレートは意図的に用意されておらず、 必要な層や演算を担うクレートを個別に use する設計です。

よく使うユーザ向けクレートは次のあたりです。

• tenferro-tensor: 実行時のテンソル値と backend trait。
• tenferro-cpu: CPU backend。
• tenferro-runtime: concrete/traced tensor helper、graph のコンパイルと実行。
• tenferro-ad: eager/traced の自動微分 API。
• tenferro-einsum: einsum 縮約。
• tenferro-linalg: 行列分解などの線形代数。
• tenferro-fft: FFT。

internal を名前に含むクレートは実装用で、ユーザ向けのAPI面ではありません。 本書では現在のAPIを暗記することが目的ではないため、 実際に使う場合は公式ドキュメントで公開クレートと import path を確認します。

メモリレイアウト

tenferro の密なテンソルは column-major（列優先）です。 つまり、左端の添字がメモリ上で最も速く変わります。

例えば、論理的な 2 x 3 の行列を

[[1, 2, 3],
 [4, 5, 6]]

と書くと、column-major の1次元バッファは次の順序になります。

[1, 4, 2, 5, 3, 6]

これは ndarray の標準的な row-major（行優先）とは違います。 ndarray やNumPyの標準的な並びからtenferroへデータを渡す場合は、 どちらの順序のバッファなのかを必ず確認します。

tenferro のテンソル生成には from_vec_col_major のように 並び順が名前に入った constructor が用意されています。 自分が渡したバッファが column-major であることを確認したうえで使います。

AI agent に実装を任せる場合も、入力データが row-major なのか column-major なのか、 各軸が何を意味するのかを明示します。 この確認を省くと、shape は合っていても、物理的には転置されたデータを計算してしまうことがあります。

使う前に確認すること

tenferro は crates.io で公開されていますが、まだ 0.x 版（pre-1.0）です。 version 間で breaking change がありうるので、 本書では現在のAPIを暗記することは目的にしません。 実際に使う場合は、公式リポジトリとドキュメントを確認してから、 小さい検証用プロジェクトで試します。

AI agent に tenferro を使ったコードを書かせる場合は、少なくとも次を確認します。

• 参照した公式ドキュメントまたはリポジトリ。
• Cargo.toml に書いたクレート名（tenferro-tensor などの個別クレート名）。
• どの名前を use したか。
• Cargo.lock に記録された version（crates.io 依存の場合）。 Git 依存を使った場合はその revision。
• CPUで実行するのか、GPUで実行するのか。
• CPUとGPUの間でデータを移す場所。
• 入出力バッファの並び順。
• 各軸の意味。
• 小さい入力で期待値を手計算できるテスト。

Cargo.lock には、実際に解決されたクレートの version が記録されます。 0.x 版のクレートは version 間でAPIが変わりうるため、 後から同じコードで再実行できるかどうかは version 固定に依存します。 そのため、演習プロジェクトでも Cargo.lock を残します。

最小例

本節の例は tenferro 0.2 で動作を確認しています。 API は 0.x の間に変わりうるので、手元で動かす場合は 公式ドキュメントで現在のシグニチャを確認してください。 以下の依存を Cargo.toml に書きます。

[dependencies]
tenferro-tensor = "0.2"
tenferro-cpu = "0.2"
tenferro-runtime = "0.2"
tenferro-einsum = "0.2"
tenferro-ad = "0.2"
tenferro-linalg = "0.2"

einsum

ij,jk->ik は行列積と同じ縮約です。 入力は column-major で渡し、結果の shape が [2, 4] になることを確認します。

use tenferro_cpu::CpuBackend;
use tenferro_einsum::TensorEinsumExt;
use tenferro_tensor::Tensor;

let a = Tensor::from_vec_col_major(vec![2, 3], vec![1.0_f64; 6])?;
let b = Tensor::from_vec_col_major(vec![3, 4], vec![1.0_f64; 12])?;
let mut backend = CpuBackend::new();

let out = [&a, &b].einsum("ij,jk->ik", &mut backend)?;
assert_eq!(out.shape(), &[2, 4]);

自動微分

AdContext で勾配を計算します。 ここでは x * x の x に関する勾配を取り、結果がスカラー（rank 0）になることを確認します。

use tenferro_ad::AdContext;
use tenferro_runtime::TracedTensor;

let ad = AdContext::builder().build()?;
let x = TracedTensor::from_vec_col_major(vec![], vec![3.0_f64])?;
let loss = (&x * &x)?;
let dx = ad.grad(&loss, &x)?;
assert_eq!(dx.rank, 0);

線形代数

tenferro-linalg は traced な行列分解を提供します。 Cholesky 分解を graph にコンパイルし、CPU backend で実行します。 linalg は CPU で BLAS/LAPACK 系の provider クレートを引き込むことがあるので、 手元で feature と外部 library の要件を確認してください。

use tenferro_cpu::CpuBackend;
use tenferro_linalg::TracedTensorLinalgExt;
use tenferro_runtime::{GraphCompiler, GraphExecutor, TracedTensor};

let a = TracedTensor::from_vec_col_major(vec![2, 2], vec![4.0_f64, 2.0, 2.0, 3.0])?;
let l = a.cholesky()?;

let mut compiler = GraphCompiler::new();
let program = compiler.compile(&l)?;

let mut executor = GraphExecutor::new(CpuBackend::new());
executor.register_extension(tenferro_linalg::register_runtime)?;
let out = executor.run(&program)?;
assert_eq!(out.shape(), &[2, 2]);

本章での位置づけ

本章では、まず標準的で軽量な ndarray を確認します。 LAPACK系の線形代数が必要なら ndarray-linalg を確認します。 そのうえで、自動微分、GPU実行、テンソル縮約のような機能が必要になった場合の 発展的な選択肢として tenferro を紹介します。

参照:

• https://github.com/tensor4all/tenferro-rs
• https://tensor4all.org/tenferro-rs/
• https://crates.io/crates/tenferro-tensor （workspace の各クレート）

データ保存とmetadata

計算結果は、後から確認できる形で保存します。画面への出力だけでは、条件を変えて比較したり、同じデータから図を作り直したり、他の言語で解析したりすることが難しくなります。

表形式データ

時系列や少数列の結果は、CSVで保存すると扱いやすいです。

t,energy,magnetization
0.0,-1.24,0.83
0.1,-1.25,0.82

CSVはgnuplot、Python、Julia、表計算ソフトなどから読みやすいため、最初の保存形式として適しています。

多次元配列

格子上の場、波動関数、画像、相関関数などの多次元配列では、CSVはあまり適していません。形状、軸、データ型を別に管理する必要があり、ファイルサイズも大きくなりやすいためです。

このような場合は、次の形式を検討します。

• HDF5
• NumPy形式（.npy, .npz）
• 必要に応じた独自binary形式

HDF5やNumPy形式は、多くの言語やツールから読み書きしやすい形式です。Rust側で計算し、保存したデータをPythonやJuliaで解析する流れにも向いています。

metadata

数値データだけでなく、実行条件も保存します。例えば、次のような情報です。

• 入力パラメータ
• 格子サイズ、shape、axisの意味
• 時間刻み、空間刻み
• 境界条件
• 乱数seed
• crate version、Rust version
• 出力ファイル名
• 実行日時

metadata は、JSON、TOML、YAMLなどのテキスト形式で保存すると確認しやすいです。

{
  "model": "ising_2d",
  "lattice_size": [64, 64],
  "temperature": 2.269,
  "seed": 12345,
  "output": "samples.npy"
}

compute と plot を分ける

本書では、計算、保存、可視化を分けます。

1. Rustで計算する。
2. 結果とmetadataを保存する。
3. 保存済みデータから図を作る。
4. 図と数値を使って検証する。

この流れにすると、AI agent に入出力コードやプロット用スクリプトを任せやすくなります。ただし、保存する列、axis、単位、境界条件、検証対象が正しいかは人間が確認します。

作業用ディレクトリのセットアップ

第2部以降の演習では、教材リポジトリの外に作業用ディレクトリを作ります。 教材リポジトリは読む場所、作業用ディレクトリは実装、テスト、結果保存を行う場所として分けます。

mkdir -p ~/rust-computational-physics-work
cd ~/rust-computational-physics-work

標準ディレクトリ構成

作業用ディレクトリの下には、章やテーマごとにsub directoryを作ります。 各leaf directoryは、原則として独立したCargo projectにします。

rust-computational-physics-work/
  AGENTS.md
  CLAUDE.md
  README.md
  calculus/
    integration/
      Cargo.toml
      src/
        lib.rs
        main.rs
      tests/
        integration_test.rs
      output/
  ode/
    euler-rk4/
  monte-carlo/
    integration/

トップレベルの指示ファイル

AI coding agentを使う場合は、作業用ディレクトリのtopにAGENTS.mdを置くことを推奨します。 Claude Codeを使う場合は、必要に応じてCLAUDE.mdも置き、AGENTS.mdを参照させます。

AGENTS.md の例:

# AGENTS.md

このディレクトリは rust-computational-physics-tutorial の演習用作業場所です。

- 教材リポジトリ本体は編集しない。
- 各テーマは独立したCargo projectとしてsub directoryに作る。
- 数値計算の本体は `src/lib.rs` に置く。
- 実行用コード、CSV出力、パラメータ実験は `src/main.rs` に置く。
- 小さい手計算可能なケースを `tests/` に追加する。
- `cargo test` を通してから結果を信用する。
- 生成されたCSVや画像は `output/` に置く。

CLAUDE.md を置く場合は、次のように AGENTS.md を参照するだけで構いません。

@AGENTS.md

テーマごとのCargo project

各テーマでは、次の形を標準にします。

mkdir -p ~/rust-computational-physics-work/calculus/integration
cd ~/rust-computational-physics-work/calculus/integration
cargo init --bin
mkdir -p tests output

• src/lib.rs: 数値計算の本体。小さい関数に分け、テストしやすくする。
• src/main.rs: 実行条件の設定、CSV出力、簡単なパラメータ走査を担当する。
• tests/: 手計算可能な入力、既知解、保存量、再現性を確認する。
• output/: CSV、metadata、図などの生成物を置く。

本文のサンプルコードも、この構成を意識して読みます。つまり、計算本体は src/lib.rs に置ける関数として定義し、fn main() はその関数を呼び出して 結果を表示したり、CSVへ保存したりする薄い実行例に留めます。

演習の標準形

第2部と第3部の演習では、各テーマに対して次の2種類を基本にします。

1. ユニットテストの追加 本文で示した関数に対して、cargo testで確認するテストを追加します。 どのケースをテストすべきかはAI coding agentと相談して構いませんが、 解析解、既知解、小さい格子、固定seed、境界条件、許容誤差を必ず確認します。
2. コードの拡張 収束表のCSV出力、保存量の記録、metadataの保存、別手法との比較などを追加します。

AI coding agentに依頼する場合も、まず関数分割、入力、出力、検証方法を確認させてから実装に進みます。

数値微分と数値積分

Important

本章の前提知識

本章の内容を理解するためには、以下の章を事前に学習しておく必要があります。

• 計算機の基本モデル
• Rustで数値計算を書く最小セット
• 多次元データと配列

第1部では、計算機の基本モデル、Rustで数値計算を行うための最小限の書き方、 多次元データ、結果保存の方針を確認しました。本章では、その知識を活用して、 微積分学の基本的な演算である「微分」と「積分」をコンピュータ上で実現する手法を学びます。

微分と積分は、物理学における運動方程式の解析、エネルギー計算、確率分布の取り扱いなど、計算物理学のあらゆる場面で必要不可欠な操作です。しかし、解析的に計算できる場合は限られており、多くの実問題では数値的な手法に頼らざるを得ません。本章では、そうした数値計算手法の原理と実装を、Rustを用いて体系的に学んでいきます。

本章の構成

• 数値微分 関数の微分をコンピュータ上で近似的に求める「数値微分」の手法を解説します。前進差分、後退差分、中心差分などの基本的な差分法と、それらの精度や誤差の特性について学びます。

• 数値積分（台形則・シンプソン則） 定積分を数値的に計算するための基本的な手法である「台形則」と「シンプソン則」を紹介します。これらは関数を区間ごとに多項式で近似し、その面積を計算することで積分値を求める古典的かつ実用的な手法です。

• ガウス求積法 台形則やシンプソン則よりも高精度な積分を実現する「ガウス求積法」について解説します。この手法は、積分点の配置を最適化することで、少ない評価点数で高い精度を達成できる強力な手法です。

• 適応型積分 被積分関数の振る舞いに応じて自動的に積分点の配置を調整する「適応型積分」を学びます。この手法により、関数が急激に変化する領域でも効率的かつ高精度に積分を計算できるようになります。

本章を終えることで、読者は微分・積分を数値的に計算する様々な手法を習得し、 それらをRustで実装する能力を身につけます。これらの技術は、後続章で扱う 微分方程式の数値解法や、より高度な物理シミュレーションの基盤となります。

作業テーマ

作業ディレクトリは ~/rust-computational-physics-work/ からの相対パスです。 ユニットテストで確認する具体的なケースは、解析解、境界条件、許容誤差の観点からAI coding agentと相談して決めます。

検証と実装の観点

本章の実装では、公式をコードに写すだけでなく、既知の関数で精度を確認します。 例えば、 の微分、低次多項式の積分、解析解が分かる滑らかな関数を使うと、 差分幅や分割数を変えたときの誤差の変化を追跡できます。

• 微分では、前進差分・中心差分の誤差が刻み幅にどう依存するかを確認する。
• 積分では、n = 0、分割数の偶奇、非滑らかな関数などの境界条件を明示する。
• differentiate、integrate、estimate_error のように小さい関数へ分ける。
• agentに実装を任せた場合も、比較対象、許容誤差、分割数をdiffで確認する。

数値微分

解析的に導関数（微分）を求めることが困難、あるいは不可能な関数に対し、計算機を用いて近似的に微分係数を求める手法を数値微分 (Numerical Differentiation) と呼びます。

物理シミュレーションにおいては、物体の運動方程式 を解く際や、ポテンシャルエネルギー から力 を求める際など、至る所で微分計算が必要となります。

本節では、数値微分の最も基本的な手法である差分法 (Finite Difference Method) について、その原理とRustによる実装、そして数値計算特有の誤差の振る舞いについて解説します。

微分の定義と差分近似

関数 の における微分係数 の数学的な定義は以下の通りです。

コンピュータでは無限小の極限 を直接扱うことはできません。その代わり、十分に小さな有限の値 （刻み幅、ステップサイズ）を用いて、微分係数を近似します。これを差分近似と呼びます。

近似の方法によって、いくつかの種類があります。

1. 前進差分 (Forward Difference)

定義式をそのまま用いて、点 とその少し先の点 の値を使う方法です。

2. 後退差分 (Backward Difference)

点 とその少し手前の点 の値を使う方法です。

3. 中心差分 (Central Difference)

点 を中心とした前後 と の値を使う方法です。

幾何学的には、前進差分と後退差分が「2点間の直線の傾き」で接線の傾きを近似するのに対し、中心差分は「前後2点の中点における傾き」を用いていることになります。一般に、中心差分の方が前進・後退差分よりも高精度になります（理由は後述します）。

Rustによる実装

これらの差分法をRustで実装してみましょう。ここでは、関数 を例にとり、 における微分係数を計算します。解析解は なので、これと比較して精度を確認します。

fn forward_difference<F>(f: F, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    (f(x + h) - f(x)) / h
}

fn backward_difference<F>(f: F, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    (f(x) - f(x - h)) / h
}

fn central_difference<F>(f: F, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    (f(x + h) - f(x - h)) / (2.0 * h)
}

fn main() {
    // 微分対象の関数 f(x) = sin(x)
    let f = |x: f64| x.sin();

    // 微分を計算する点
    let x: f64 = 1.0;

    // 刻み幅 h
    let h: f64 = 0.01;

    // 解析解 (真値): cos(1.0)
    let exact = x.cos();

    // 1. 前進差分
    let diff_forward = forward_difference(f, x, h);

    // 2. 後退差分
    let diff_backward = backward_difference(f, x, h);

    // 3. 中心差分
    let diff_central = central_difference(f, x, h);

    println!("解析解:    {:.10}", exact);
    println!("前進差分:  {:.10} (誤差: {:.2e})", diff_forward, (diff_forward - exact).abs());
    println!("後退差分:  {:.10} (誤差: {:.2e})", diff_backward, (diff_backward - exact).abs());
    println!("中心差分:  {:.10} (誤差: {:.2e})", diff_central, (diff_central - exact).abs());
}

実行結果の例：

解析解:    0.5403023059
前進差分:  0.5360859810 (誤差: 4.22e-3)
後退差分:  0.5445006207 (誤差: 4.20e-3)
中心差分:  0.5402933009 (誤差: 9.00e-6)

この結果から、同じ刻み幅 でも、中心差分の方が圧倒的に誤差が小さいことがわかります。

誤差の解析：なぜ中心差分が良いのか？

テイラー展開を用いて誤差（打ち切り誤差）を評価することで、この理由を数学的に説明できます。

関数 を の周りでテイラー展開すると以下のようになります。

前進差分の誤差

前進差分の式にテイラー展開を代入して整理すると：

つまり、近似値と真値 との差（誤差）は、 の1乗に比例します（オーダー ）。

中心差分の誤差

一方、中心差分の式 に代入すると、偶数次の項が相殺されます。

こちらの誤差は、 の2乗に比例します（オーダー ）。

が小さい場合（例： ）、 となり、 に比べてはるかに小さくなります。これが、中心差分が高精度である理由です。

最適なステップサイズと丸め誤差

「 を小さくすればするほど精度が良くなる」と考えがちですが、数値計算ではそう単純ではありません。コンピュータの浮動小数点数には丸め誤差が存在するためです。

を極端に小さくすると、分子の において、非常に近い値同士の引き算（桁落ち）が発生し、有効桁数が失われてしまいます。

実際に を変化させながら誤差をプロットしてみると、以下のようになります。

1. 打ち切り誤差支配領域: が大きい領域。テイラー展開の近似誤差が支配的で、 を小さくすると精度が向上する。
2. 丸め誤差支配領域: が小さすぎる領域。桁落ちによる丸め誤差が支配的で、 を小さくすると逆に精度が悪化する。

最適な は、これら2つの誤差のバランスが取れる点に存在します。

実験：最適なhを探る

以下のコードで、 を から まで変化させたときの誤差を確認してみましょう。

fn central_difference<F>(f: F, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    (f(x + h) - f(x - h)) / (2.0 * h)
}

fn central_difference_error<F, G>(f: F, df_exact: G, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
    G: Fn(f64) -> f64,
{
    (central_difference(f, x, h) - df_exact(x)).abs()
}

fn main() {
    let f = |x: f64| x.sin();
    let df_exact = |x: f64| x.cos();
    let x: f64 = 1.0;

    println!("h,      Error (Central)");

    let mut h = 1.0;
    for _ in 0..16 {
        h /= 10.0;
        let error = central_difference_error(f, df_exact, x, h);
        println!("{:.1e}, {:.2e}", h, error);
    }
}

実行結果：

h,      Error (Central)
1.0e-1, 9.00e-4
1.0e-2, 9.00e-6
1.0e-3, 9.01e-8
1.0e-4, 9.00e-10
1.0e-5, 1.11e-11
1.0e-6, 2.77e-11
1.0e-7, 1.94e-10
1.0e-8, 2.97e-9
...
1.0e-16, 1.48e-2

倍精度浮動小数点数（f64）を用いる場合、中心差分における最適な はおおよそ 〜 程度になります（関数のスケールにもよりますが、マシンイプシロン に対して 程度が目安と言われます）。

Note

実用的な数値計算ライブラリでは、このような最適な を自動的に推定するアルゴリズムが組み込まれている場合もあります。

2階微分の数値計算

運動方程式 などで登場する2階微分 も、中心差分法を応用して求めることができます。

この近似式の誤差も です。

fn second_derivative_central<F>(f: F, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    (f(x + h) - 2.0 * f(x) + f(x - h)) / (h * h)
}

まとめ

• 数値微分には前進差分、後退差分、中心差分などの手法がある。
• 中心差分は の精度を持ち、他の2つ（ ）より高精度であるため、通常はこれが推奨される。
• ステップサイズ は小さければ良いわけではなく、丸め誤差（桁落ち） に注意する必要がある。
• f64 を用いる場合、中心差分での最適な は 程度が目安となる。

次節では、数値積分について学びます。

数値積分

数値積分（Numerical Integration）は、定積分

の値を数値的に近似計算する手法です。解析的に不定積分を求めることが困難な関数や、離散的なデータ点としてしか与えられていない関数（実験データなど）の積分値を求めるために用いられます。物理学では、運動方程式の積分によるエネルギー計算、電磁場のポテンシャル計算、分配関数の計算など、極めて広範な応用を持ちます。

本節では、最も基本的かつ実用的な手法である台形則 (Trapezoidal Rule) と シンプソン則 (Simpson’s Rule) について解説します。これらの手法は「ニュートン・コーツ (Newton-Cotes) の公式」と呼ばれる手法群の一種で、積分区間を等間隔に分割して関数値を評価します。

複合台形則 (Trapezoidal Rule)

原理

台形則は、関数 を各小区間で一次関数（直線）で近似し、その下の面積を台形の面積として計算する手法です。

積分区間 を 等分し、刻み幅を とします。分点を ( ) とすると、各区間 の面積は以下のように近似されます。

これらを全区間にわたって足し合わせることで、全体の積分値の近似式（複合台形則）が得られます。

誤差

台形則の誤差は、刻み幅 の2乗に比例します（ ）。つまり、分割数 を2倍（ を半分）にすると、誤差は約 1/4 に減少します。

複合シンプソン則 (Simpson’s Rule)

原理

シンプソン則は、隣り合う2つの小区間 をまとめて、その区間内の関数を二次関数（放物線）で近似する手法です。これにより、一次近似である台形則よりも高い精度が得られます。

シンプソン則を適用するには、分割数 が偶数である必要があります。近似式は以下のようになります。

係数が のパターンで並ぶのが特徴です。

誤差

シンプソン則の誤差は、刻み幅 の4乗に比例します（ ）。つまり、分割数 を2倍にすると、誤差は約 1/16 に激減します。これが、シンプソン則が広く使われる理由です。

Rustによる実装

それでは、これらの手法をRustで実装し、精度を比較してみましょう。 例として、 を区間 で積分します。解析解は です。

use std::f64::consts::PI;

/// 台形則
fn trapezoidal_rule<F>(f: F, a: f64, b: f64, n: usize) -> f64
where
    F: Fn(f64) -> f64,
{
    let h = (b - a) / n as f64;
    let sum: f64 = (1..n).map(|i| f(a + i as f64 * h)).sum();

    // 両端の点は重み 1/2
    h * (0.5 * f(a) + sum + 0.5 * f(b))
}

/// シンプソン則 (nは偶数でなければならない)
fn simpsons_rule<F>(f: F, a: f64, b: f64, n: usize) -> f64
where
    F: Fn(f64) -> f64,
{
    assert!(n % 2 == 0, "Simpson's rule requires an even number of intervals.");

    let h = (b - a) / n as f64;

    let mut sum_odd = 0.0;
    let mut sum_even = 0.0;

    for i in 1..n {
        let x = a + i as f64 * h;
        if i % 2 == 0 {
            sum_even += f(x);
        } else {
            sum_odd += f(x);
        }
    }

    h / 3.0 * (f(a) + 4.0 * sum_odd + 2.0 * sum_even + f(b))
}

fn main() {
    let f = |x: f64| x.sin();
    let a = 0.0;
    let b = PI;
    let exact = 2.0;

    println!("{:<10} {:<20} {:<20} {:<20} {:<20}",
        "N", "Trapezoidal", "Error (Trap)", "Simpson", "Error (Simp)");
    println!("{}", "-".repeat(95));

    let n_values = [10, 20, 40, 80, 160];

    for &n in &n_values {
        let trap = trapezoidal_rule(f, a, b, n);
        let simp = simpsons_rule(f, a, b, n);

        println!(
            "{:<10} {:.15}    {:.5e}           {:.15}    {:.5e}",
            n, trap, (trap - exact).abs(), simp, (simp - exact).abs());
    }
}

実行結果：

N          Trapezoidal          Error (Trap)         Simpson              Error (Simp)
-----------------------------------------------------------------------------------------------
10         1.983523537509455    1.64765e-2           2.000109517315004    1.09517e-4
20         1.995885972708715    4.11403e-3           2.000006784441801    6.78444e-6
40         1.998971810497066    1.02819e-3           2.000000423093183    4.23093e-7
80         1.999742972445836    2.57028e-4           2.000000026428759    2.64288e-8
160        1.999935744350136    6.42556e-5           2.000000001651570    1.65157e-9

結果の考察

• 台形則: が2倍（10 -> 20）になると、誤差は約 1/4 ( ) になっています。これは理論通り の収束です。
• シンプソン則: が2倍（10 -> 20）になると、誤差は約 1/16 ( ) になっています。これも理論通り の収束を示しています。
• 同じ分割数 で比較すると、シンプソン則の方が圧倒的に高精度であることがわかります。

どちらを使うべきか？

• 関数が滑らかな場合: シンプソン則が推奨されます。少ない計算コストで高い精度が得られます。
• 実験データなどノイズを含む場合: 高次の近似がかえってノイズを増幅させる可能性があるため、台形則の方が安定する場合があり、一般的によく用いられます。
• 周期関数の全周期積分: 周期関数の場合、台形則が例外的に極めて高い精度（指数関数的収束）を示すことが知られています。

まとめ

• 台形則は一次近似であり、誤差は 。シンプルでロバスト。
• シンプソン則は二次近似であり、誤差は 。滑らかな関数に対して非常に高精度。
• シンプソン則を使うには、分割数 が偶数である必要がある。

次節では、より高度な積分手法である「ガウス求積法」について学びます。

ガウス求積法

前節で紹介した台形則やシンプソン則は、積分区間を等間隔に分割する手法（ニュートン・コーツの公式）でした。これに対し、ガウス求積法 (Gaussian Quadrature)、特にガウス・ルジャンドル積分 (Gauss-Legendre Integration) は、積分点（分点）の配置と重みを最適化することで、より少ない点数で極めて高い精度を実現する手法です。

この手法は、被積分関数が滑らかである（高階微分可能である）場合に圧倒的な威力を発揮します。

原理

基本公式

ガウス・ルジャンドル積分では、積分区間を に正規化した上で、以下の近似式を用います。

ここで、 は積分点の数（次数）、 は積分点、 は重み係数です。 ガウス求積法の最大の特徴は、 個の点を用いるだけで、 次以下の任意の多項式を厳密に積分できる という点です。

例えば、 の場合、5次関数までの積分が厳密に一致します。シンプソン則（3点等間隔）が3次関数までしか厳密でないことと比較すると、その効率の良さが分かります。

積分点と重み

積分点 は、ルジャンドル多項式 (Legendre Polynomials) の零点（根）として選ばれます。 また、重み は以下の式で決定されます。

これらは解析的に計算可能であり、多くの数値計算ライブラリではあらかじめ計算された値（テーブル）が用いられます。

例： n=5 の場合の積分点と重み

積分点は原点に対して対称に配置され、端点（ ）よりも内側に配置されるのが特徴です。

積分区間の変換

一般的な積分区間 での積分計算を行うためには、変数変換を行って積分区間を に写像する必要があります。

変数変換：

これにより、公式は以下のように変形されます。