matlab-scientific-computing

MATLAB/Octave — 科学技術計算

概要

MATLAB は、行列演算と科学技術計算に最適化された数値計算環境です。GNU Octave は、互換性の高い無料のオープンソース代替ソフトウェアです。すべてのコード例では、両方のプラットフォームで実行できる MATLAB 構文を使用しています。

使用する場面

• 行列演算や線形代数（固有値、SVD、最小二乗法）の実行
• 常微分方程式や偏微分方程式の数値解法
• 信号処理（FFT、フィルタリング、スペクトル解析）
• 2D/3D の科学技術向け可視化や論文用の図の作成
• 数値最適化や求根
• 統計解析や曲線フィッティング
• 実験データファイルのバッチ処理
• Python ベースの数値計算には、代わりに numpy/scipy を使用してください
• 推論を伴う統計モデリングには、代わりに statsmodels を使用してください

前提条件

# GNU Octave (free, open-source)
# macOS
brew install octave
# Ubuntu/Debian
sudo apt install octave

# Running scripts
octave script.m                           # Octave
matlab -nodisplay -nosplash -r "run('script.m'); exit;"  # MATLAB

注記：MATLAB の使用には、MathWorks の商用ライセンスが必要です。GNU Octave は無料で、ほとんどの MATLAB スクリプトを変更せずに実行できます。Octave との主な違い：# コメントと ++／+= 演算子をサポートしています。一部の MATLAB ツールボックス関数は使用できません。

クイックスタート

% Load data, fit, and plot
x = linspace(0, 2*pi, 100);
y = sin(x) + 0.1 * randn(size(x));
p = polyfit(x, y, 5);
y_fit = polyval(p, x);

figure;
plot(x, y, 'bo', x, y_fit, 'r-', 'LineWidth', 2);
xlabel('x'); ylabel('y');
legend('Data', 'Polynomial fit');
title('Curve Fitting Example');
saveas(gcf, 'fit_result.png');

コア API

1. 行列演算

MATLAB は、基本的に行列と配列を扱います。

% Create matrices
A = [1 2 3; 4 5 6; 7 8 9];    % 3x3 matrix
v = linspace(0, 1, 100);       % 100 evenly spaced points
I = eye(3);                     % Identity matrix
R = rand(3, 3);                 % Uniform random
N = randn(3, 3);                % Normal random

% Operations
B = A';                  % Transpose
C = A * B;               % Matrix multiplication
D = A .* B;              % Element-wise multiplication
x = A \ [1; 2; 3];      % Solve Ax = b (preferred over inv(A)*b)
fprintf('Solution: [%.2f, %.2f, %.2f]\n', x);

% Indexing and manipulation
A = magic(5);
sub = A(1:3, 2:4);      % Submatrix (rows 1-3, cols 2-4)
row = A(2, :);           % Entire row 2
col = A(:, 3);           % Entire column 3
A(A < 5) = 0;           % Logical indexing

% Concatenation
C = [A; ones(1, 5)];    % Vertical (add row)
D = [A, zeros(5, 1)];   % Horizontal (add column)
fprintf('Size: %d x %d\n', size(C));

2. 線形代数

A = [4 1 2; 1 3 1; 2 1 5];

% Eigendecomposition
[V, D] = eig(A);           % V: eigenvectors, D: diagonal eigenvalues
fprintf('Eigenvalues: %.2f, %.2f, %.2f\n', diag(D));

% Singular value decomposition
[U, S, V] = svd(A);
fprintf('Singular values: %.2f, %.2f, %.2f\n', diag(S));

% Matrix decompositions
[L, U, P] = lu(A);         % LU with pivoting
[Q, R] = qr(A);            % QR decomposition
R_chol = chol(A);           % Cholesky (symmetric positive definite)

% Condition number and rank
fprintf('Condition number: %.2f\n', cond(A));
fprintf('Rank: %d\n', rank(A));

3. プロットと可視化

% 2D line plots
x = 0:0.1:2*pi;
figure;
plot(x, sin(x), 'b-', 'LineWidth', 2); hold on;
plot(x, cos(x), 'r--', 'LineWidth', 2);
xlabel('x'); ylabel('y');
title('Trigonometric Functions');
legend('sin(x)', 'cos(x)');
grid on;
saveas(gcf, 'trig.png');

% 3D surface plot
[X, Y] = meshgrid(-2:0.1:2, -2:0.1:2);
Z = X.^2 + Y.^2;
figure;
surf(X, Y, Z);
colorbar; xlabel('X'); ylabel('Y'); zlabel('Z');
title('Paraboloid');
print('-dpdf', 'surface.pdf');

% Multi-panel figure
figure;
subplot(2, 2, 1); plot(x, sin(x)); title('sin');
subplot(2, 2, 2); plot(x, cos(x)); title('cos');
subplot(2, 2, 3); bar([1 3 2 5 4]); title('Bar');
subplot(2, 2, 4); histogram(randn(1000, 1), 30); title('Histogram');
saveas(gcf, 'panels.png');

4. データのインポート／エクスポート

% CSV / tabular data
T = readtable('data.csv');
M = readmatrix('data.csv');
fprintf('Table: %d rows x %d cols\n', height(T), width(T));

% Write data
writetable(T, 'output.csv');
writematrix(M, 'output.csv');

% MAT files (MATLAB native binary)
A = rand(100, 100);
save('data.mat', 'A');          % Save variable
S = load('data.mat', 'A');     % Load specific variable

% Images
img = imread('image.png');
fprintf('Image size: %d x %d x %d\n', size(img));
imwrite(img, 'output.jpg');

5. 統計とデータ分析

data = randn(1000, 1) * 5 + 50;

% Descriptive statistics
fprintf('Mean: %.2f, Std: %.2f, Median: %.2f\n', mean(data), std(data), median(data));
fprintf('Min: %.2f, Max: %.2f\n', min(data), max(data));

% Correlation and covariance
X = randn(100, 3);
R = corrcoef(X);
fprintf('Correlation matrix:\n');
disp(R);

% Linear regression (polyfit)
x = (1:50)';
y = 2.5 * x + 10 + randn(50, 1) * 5;
p = polyfit(x, y, 1);
fprintf('Slope: %.2f, Intercept: %.2f\n', p(1), p(2));

% Moving statistics
y_smooth = movmean(y, 5);

6. 微分方程式

% First-order ODE: dy/dt = -2y, y(0) = 1
f = @(t, y) -2 * y;
[t, y] = ode45(f, [0 5], 1);
figure; plot(t, y, 'b-', 'LineWidth', 2);
xlabel('Time'); ylabel('y(t)');
title('Exponential Decay');
fprintf('Final value: %.4f (expected: %.4f)\n', y(end), exp(-10));

% Second-order ODE: y'' + 0.5y' + 4y = 0 (damped oscillator)
% Convert to system: y1' = y2, y2' = -0.5*y2 - 4*y1
f = @(t, y) [y(2); -0.5*y(2) - 4*y(1)];
[t, y] = ode45(f, [0 20], [1; 0]);
figure; plot(t, y(:,1), 'b-', 'LineWidth', 2);
xlabel('Time'); ylabel('Displacement');
title('Damped Oscillator');

7. 信号処理

% Generate signal with two frequencies
fs = 1000;                       % Sampling frequency
t = 0:1/fs:1-1/fs;
signal = sin(2*pi*50*t) + 0.5*sin(2*pi*120*t) + randn(size(t))*0.2;

% FFT
Y = fft(signal);
f = (0:length(Y)-1) * fs / length(Y);
figure;
plot(f(1:length(f)/2), abs(Y(1:length(Y)/2)));
xlabel('Frequency (Hz)'); ylabel('|FFT|');
title('Frequency Spectrum');

% FIR low-pass filter (keep < 80 Hz)
b = fir1(50, 80/(fs/2));          % 50th order, cutoff 80 Hz
filtered = filter(b, 1, signal);
figure;
plot(t, signal, 'b', t, filtered, 'r', 'LineWidth', 1.5);
legend('Original', 'Filtered');
title('Low-pass Filtering');

8. 関数とプログラミング

% Anonymous functions
f = @(x) x.^2 + 2*x + 1;
fprintf('f(5) = %d\n', f(5));   % 36

% Function files (save as myfunc.m)
% function [rmse, r2] = myfunc(y_true, y_pred)
%     residuals = y_true - y_pred;
%     rmse = sqrt(mean(residuals.^2));
%     ss_res = sum(residuals.^2);
%     ss_tot = sum((y_true - mean(y_true)).^2);
%     r2 = 1 - ss_res / ss_tot;
% end

% Struct for organized data
experiment.name = 'Trial 1';
experiment.data = rand(10, 3);
experiment.params = struct('temp', 37, 'pH', 7.4);
fprintf('Experiment: %s, %d samples\n', experiment.name, size(experiment.data, 1));

主要な概念

ベクトル化

MATLAB はベクトル化演算に最適化されています。可能な場合は、明示的なループを避けてください。

% Slow (loop)
n = 1e6;
y = zeros(1, n);
tic;
for i = 1:n
    y(i) = sin(i/1000);
end
t_loop = toc;

% Fast (vectorized)
tic;
y = sin((1:n)/1000);
t_vec = toc;
fprintf('Loop: %.3fs, Vectorized: %.3fs, Speedup: %.1fx\n', t_loop, t_vec, t_loop/t_vec);

要素ごとの演算と行列演算

演算子	行列	要素ごと
乗算	A * B	A .* B
除算	A / B（右除算）	A ./ B
べき乗	A ^ n（行列のべき乗）	A .^ n

ODE ソルバーの選択

ソルバー	次数	使用する場面
ode45	4-5	デフォルト。ほとんどの非スティッフ問題
ode23	2-3	精度が粗い解を求める場合。1ステップあたりの計算が高速
ode113	可変	高精度が必要で、評価コストが高い場合
ode15s	可変	スティッフ問題（化学反応速度論、回路）
ode23s	2	スティッフ問題で、中程度の精度が必要な場合

よくあるワークフロー

ワークフロー：データ分析パイプライン

% 1. Load data
data = readtable('experiment.csv');

% 2. Clean
data = rmmissing(data);
fprintf('After cleaning: %d rows\n', height(data));

% 3. Group analysis
groups = unique(data.Category);
results = table();
for i = 1:length(groups)
    mask = strcmp(data.Category, groups{i});
    subset = data(mask, :);
    row = table(groups(i), mean(subset.Value), std(subset.Value), ...
        'VariableNames', {'Category', 'Mean', 'Std'});
    results = [results; row];
end
disp(results);

% 4. Visualize and save
figure;
bar(categorical(results.Category), results.Mean);
hold on;
errorbar(1:height(results), results.Mean, results.Std, '.k');
ylabel('Mean Value'); title('Results by Category');
saveas(gcf, 'results.png');
writetable(results, 'summary.csv');

ワークフロー：数値シミュレーション（熱方程式）

% 1D heat equation: du/dt = alpha * d2u/dx2
L = 1; N = 100; T_end = 0.5; alpha = 0.01;
dx = L / (N - 1);
dt = 0.4 * dx^2 / alpha;  % CFL condition
x = linspace(0, L, N);
nsteps = floor(T_end / dt);

% Initial condition: Gaussian pulse
u = exp(-50 * (x - 0.5).^2);

% Time stepping (explicit finite difference)
figure;
for step = 1:nsteps
    u_new = u;
    for i = 2:N-1
        u_new(i) = u(i) + alpha * dt / dx^2 * (u(i+1) - 2*u(i) + u(i-1));
    end
    u = u_new;
    if mod(step, floor(nsteps/5)) == 0
        plot(x, u, 'LineWidth', 1.5); hold on;
    end
end
xlabel('Position'); ylabel('Temperature');
title('Heat Equation Evolution');
legend(arrayfun(@(n) sprintf('t=%.2f', n*dt*floor(nsteps/5)), 1:5, 'UniformOutput', false));
saveas(gcf, 'heat_equation.png');

ワークフロー：ファイルのバッチ処理

1. dir('data/*.csv') でファイルを一覧表示します
2. ファイルを順に処理し、それぞれを readtable() で読み込みます
3. 各ファイルに分析関数を適用します
4. vertcat() で結果を要約テーブルにまとめます
5. writetable() で要約をエクスポートします

主要なパラメーター

パラメーター	関数	デフォルト	範囲	効果
次数	polyfit	—	1–20	フィッティングに使用する多項式の次数
tspan	ode45	—	[t0, tf]	積分する時間区間
'LineWidth'	plot	0.5	0.1–5	プロットの線の太さ
bins	histogram	自動	1–1000	ヒストグラムのビン数
フィルター次数	fir1	—	10–200	FIR フィルターの次数（大きいほど急峻）
'RelTol'	ode45	1e-3	1e-12–1e-1	相対誤差の許容値
'AbsTol'	ode45	1e-6	1e-15–1e-1	絶対誤差の許容値

ベストプラクティス

1. ループより常にベクトル化を優先する：MATLAB の JIT も優れていますが、大きな配列ではベクトル化したコードの方が 10〜100 倍高速です。bsxfun、論理インデックス、配列演算を使用してください。

2. 配列を事前に確保する：ループ内で配列を拡張すると、メモリの割り当てが繰り返し発生します。

   % Bad: y = []; for i=1:n, y = [y, f(i)]; end
   % Good:
   y = zeros(1, n);
   for i = 1:n, y(i) = f(i); end
   

3. 連立一次方程式には inv() ではなく \ を使用する：A\b は inv(A)*b より数値的に安定しており、高速です。

4. アンチパターン — 浮動小数点数の比較に == を使用する：浮動小数点数の精度を考慮し、a == b ではなく abs(a - b) < tol を使用してください。

5. 論文用の図はベクター形式で保存する：拡大縮小可能な図には、PNG ではなく print('-dpdf', 'fig.pdf') または print('-dsvg', 'fig.svg') を使用してください。

6. アンチパターン — 0 始まりと 1 始まりのロジックを混在させる：MATLAB の配列は 1 から始まります。Python/C から移植する場合は、すべてのインデックスを調整してください。

7. 書式付き出力には disp より fprintf を使用する：fprintf なら数値の書式を制御できますが、disp は値をそのまま表示するだけです。

よく使うレシピ

レシピ：信頼区間付きの曲線フィッティング

x = (1:20)';
y = 3 * exp(-0.2 * x) + 0.5 * randn(size(x));

% Nonlinear fit using lsqcurvefit (Optimization Toolbox)
% Or use polyfit for polynomial:
p = polyfit(x, y, 3);
y_fit = polyval(p, x);
residuals = y - y_fit;
rmse = sqrt(mean(residuals.^2));
fprintf('RMSE: %.4f\n', rmse);

figure;
plot(x, y, 'ko', x, y_fit, 'r-', 'LineWidth', 2);
xlabel('x'); ylabel('y');
title(sprintf('Polynomial Fit (RMSE = %.3f)', rmse));
saveas(gcf, 'curvefit.png');

レシピ：画像処理の基本

img = imread('sample.png');
gray = rgb2gray(img);

% Edge detection
edges = edge(gray, 'Canny');

% Display
figure;
subplot(1, 3, 1); imshow(img); title('Original');
subplot(1, 3, 2); imshow(gray); title('Grayscale');
subplot(1, 3, 3); imshow(edges); title('Edges');
saveas(gcf, 'image_processing.png');

レシピ：Python との連携

% Call Python from MATLAB (requires Python on PATH)
result = py.numpy.array([1, 2, 3, 4, 5]);
py_mean = py.numpy.mean(result);
fprintf('Python numpy mean: %.1f\n', double(py_mean));

% Call MATLAB from Python (requires matlab.engine)
% import matlab.engine
% eng = matlab.engine.start_matlab()
% result = eng.sqrt(42.0)

トラブルシューティング

問題	原因	解決策
Undefined function or variable	関数がパス上にない、または名前の誤り	which funcname で確認し、addpath() でパスを追加します
* の次元が一致しない	行列のサイズに互換性がない	要素ごとの演算には .* を使用し、両方のオペランドの size() を確認します
ODE ソルバーが非常に遅い	非スティッフ系向けのソルバーでスティッフ問題を解いている	スティッフ系には ode15s または ode23s に切り替えます
Singular matrix 警告	行列のランクが不足している	cond(A) を確認し、pinv()（擬似逆行列）を使用するか、正則化します
Octave の pkg エラー	パッケージがインストールされていない	pkg install -forge package_name; pkg load package_name を実行します
図が保存されない	有効な figure ハンドルがない	プロットの前に figure; を実行して、明示的に図を作成します
Out of memory	大きな配列の割り当て	スパース行列（sparse()）を使用する、チャンクごとに処理する、またはスワップ領域を増やします

関連スキル

• statsmodels-statistical-modeling — 推論テーブルを含む統計モデリング向けの Python の代替手段
• matplotlib-scientific-plotting — Python でのプロット。Python エコシステムで作業する場合に使用します
• scikit-learn-machine-learning — Python の ML。Python に適した分類／クラスタリングのタスク向けです

参考資料

• GNU Octave のドキュメント — 公式マニュアルと関数リファレンス
• MATLAB のドキュメント — MathWorks の公式リファレンス
• MATLAB–Octave の互換性 — 主な構文と機能の違い