目的関数のスケーリング
目的関数は１の大きさになるようにするとよい。
𝜌
𝑀 𝑑𝑆 = 1
パラメタ推定
M. Hashiguchi and D. Mi 2018 90
最小二乗法によるパラメタ推定
実験式のパラメタ推定
非線形材料の伸び
λ
と印加圧力P
の関係式が与えられているが、係数
C10
、C01
を推定したい。M. Hashiguchi and D. Mi 2018 92
モデルビルダ
𝜆 → 𝑃(𝜆)
の関係を定義する実験値を用意する
実験値の第１列に相当するパラメタが、
COMSOL
側の変数lambda
に対応すると決める
実験値の第２列が、
COMSOL
側の変数 Ｐに対応すると決める（実験値を外部ファイルから読み込む 場合は、後ほど説明する）
実験値をマニュアルで入力する。
M. Hashiguchi and D. Mi 2018 94
ソルバー
Levenberg-Marquardt
最小二乗法用のソルバー
計算結果の検証
実験値を表示する 変数を利用できる。
M. Hashiguchi and D. Mi 2018 96
伝熱モデルで材料を推定してみる
𝑇
以外の境界条件を固定し、𝑇
を変化させたときの温度を算出し、それらを実験データとする。
次に、それらのデータを使って、熱伝導係数を推定してみる。
数値実験データの作成
（１）データ採取点を下図のように決める。（ジオメトリでポイント作成）
（２） スタディ１にパラメトリックスイープを追加し、
温度
T1
を変化させる。（３） 結果：計算値で上述の
3
点を選び、ポイント評価をし、クリップボードにコピー後、エクセルにペーストする。
（４） エクセルファイルを保存しておく。
ここから
csv
ファイルを作成し、それをファイルとして読み込む。M. Hashiguchi and D. Mi 2018 98
パラメタ推定
測定点の情報を取り出す
あとでプロットする ときの名前(任意）
結果
真値
400 W/(m
・K)
に対して405.26 (1.3%
の誤差）＜面白いこと＞
外部ファイルからではなく、結果テーブルを 利用すると、
400
と正確値が算出される。従って、
Verification
を最も正確に行う場合 には、結果テーブルを利用する方が良い。テーブルのフル精度表示による出力も可。
M. Hashiguchi and D. Mi 2018 100
このことを避けるためには
次ページのフル精度表示を行った
上で、外部ファイルに吐き出せば良い。
テーブルのフル精度表示
（１） 外部ファイル（ｃｓｖ）へのテーブル出力の前に、
フル精度ボタンをクリックすることで、表示桁をフル精度に変更できる。
(2)
正確値を推定できる。セミナでの質問に関する回答
時系列データが与えられた場合
M. Hashiguchi and D. Mi 2018 102
伝熱への応用
Ed Fontes
https://www.comsol.jp/video/using-the-optimization -module-for-parameter-estimation
温度測定データ
T(t)
温度の時系列データから熱伝導係数を推定するモデルビルダ
M. Hashiguchi and D. Mi 2018 104
結果
熱伝導係数
化学工学への応用
温度の異なる
5
種類の時系列 データを読み込むM. Hashiguchi and D. Mi 2018 106
結果
計算時間
19
秒𝐴 𝐸
アプリ
M. Hashiguchi and D. Mi 2018 108
リチウム - イオン電池への応用
モデルビルダ
外部ファイル
(.csv)
を読み込むM. Hashiguchi and D. Mi 2018 110
アプリ画面
音振動とプログラミング
M. Hashiguchi and D. Mi 2018 112
アプリケーションライブラリ
ファイル～アプリケーションライブラリ
アプリによる概要説明
モデルビルダーでアプリケーションビルダーを クリックすればそちらへ移行する
アプリの実行のみを行う場合
問題説明、理論、参考文献、操作手順、結果の説明など
M. Hashiguchi and D. Mi 2018 114
すぐに計算に取り掛かれる
ソフトウェアの使い方を知らなくても良い！ 誰もが即戦力になる。
モデルの内容を見る
アプリケーションでない項目の例題は「開く」でモデルビルダーを表示する。
モデルビルダへ
モデルビルダーでアプリケーションビルダーを クリックすればそちらへ移行する
アプリの実行のみを行う場合
問題説明、理論、参考文献、操作手順、結果の説明など
M. Hashiguchi and D. Mi 2018 116
モデルビルダーの内容
最適化がない！
何故？
割線法（セカント法）
𝑓 𝑥 = 0
となる𝑥
を求める。𝑥 = 𝑥 − 𝑓(𝑥 ) 𝑥 − 𝑥
𝑓 𝑥 − 𝑓(𝑥 )
𝑓 𝑥 = 𝑓 𝑥 + 𝑥 − 𝑥 ~𝑓 𝑥 + 𝑓 (𝑥 )(𝑥 − 𝑥 )
いま、𝑓 𝑥 = 0
となったとすると𝑓 𝑥 + 𝑓 𝑥 𝑥 − 𝑥 = 0 𝑥 = 𝑥 − 𝑓(𝑥 )
𝑓 (𝑥 )
これはニュートン法であり、
関数値と関数の微分値の 両方の計算が必要である。
𝑓 𝑥 = 𝑓 𝑥 − 𝑓(𝑥 ) 𝑥 − 𝑥
後退差分で近似
𝑓
の計算のみで済む（計算速度が速い）
M. Hashiguchi and D. Mi 2018 118
アプリケーションビルダへ
ここを分析して みる
クリック
FEM を介した代数方程式
ＦＥＭの役目
構造物の腕の長さ
L
を決めて、その固有振動数ｆを計算する 関数化frequency(
引数 Ｌ)
こうすれば、長さ
L
の構造物の固有振動数frequency(L)
を目標値targetfq
に 一致させたいという要求は、見かけ上、次の代数方程式で表現できる。ただし、
𝑓(𝐿) = frequency L − targetfq 𝑓(𝐿) = 0
数値解法として、セカント法を利用すれば、Ｌを求めることができる。
COMSOL
ではアプリケーションビルダのメソッドエディタでプログラムを実装できる。
M. Hashiguchi and D. Mi 2018 120
FEM の関数化
// Function used in compute_and_play method // Syntax: double f = frequency(double L)
model.param().set("Lp", pronglength);
model.study("std1").run();
set_results(); // in case there was no solution included in the application file double[][] d = model.result().numerical("gev1").getReal();
f = d[0][0];
メソッド
frequency
Ｌｐに
pronglength
の数値をセットするFEM
の実行計算された固有周波数を
f
とするセカント法
long t0 = timeStamp(); //
計算時間の記録の初期化solution_state = "nosolution"; //
showProgress();
if (findlength) { // findlength declared variable linked to check box in main form setProgress(0, "Computing prong length.");
//
セカント法の開始:
int MAXITERATIONS = 20; //
最大反復回数double L1 = 85; //
長さの初期値20 Hz<fq<10,000 Hz.
double L2 = 60; //
長さの初期値20 Hz<fq<10,000 Hz.
double carry = L1; carry
にL1
をコピーしておくdouble f2 = frequency(L2)-targetfq; // L2
に対するＦＥＭ計算結果でf(L2)
を計算setProgress(100/MAXITERATIONS); //
プログレスバーの設定fq = frequency(L1); //
画面表示用の変数、L1
に対するFEM
計算結果をfq
とするsetProgress(20);
メソッド
b_solved_and_update_results
M. Hashiguchi and D. Mi 2018 122
セカント法（続き）
double f1 = fq-targetfq; //f(L1)
を計算L1 = L1-f1*((L1-L2)/(f1-f2)); //
セカント法L2 = carry; //caryy
の内容をL2
にコピーL1 = Math.max(L1, 1e-3); //Java
の数学ライブラリの最大値計算関数int k = 2; //
整数k
を2
とするwhile (k < MAXITERATIONS && Math.abs(f1) > fqtol) { f2 = f1; //f1
の内容をf2
に置き換えるfq = frequency(L1); //fq
をL1
でのFEM
計算結果とするf1 = fq-targetfq; //f1
を更新するcarry = L1; // carry
にL1
をコピーするL1 = L1-f1*((L1-L2)/(f1-f2)); //
新しいL1
を算出するL2 = carry; //carry
をL2
にコピーするL1 = Math.max(L1, 1e-3); //
k = k+1;
setProgress(k*100/MAXITERATIONS);
}
𝑥 = 𝑥 − 𝑓(𝑥 ) 𝑥 − 𝑥 𝑓 𝑥 − 𝑓(𝑥 )
𝐿 𝐿
𝐿 𝐿
carry secant
𝑓 𝑓
𝑓 𝑓
FEM FEM
FEM
反復計算
セカント法（続き）
L1 = Math.round(L1*100)/100.00; // we won't get more than 2 decimal accuracy (mesh limits the accuracy rather than the secant method)
model.param().set(“Lp”, L1);
グローバル定義：パラメタLp
への収束値L1
のセットfq = fq;
収束値L1
に対応する固有周波数のセットif (Math.abs(f1) > fqtol) {
error("Computation terminated after "+toString(MAXITERATIONS)+" iterations.
Frequency diff: "+toString(Math.abs(f1)));
} }
else { // "Find prong length" check box is cleared in main form setProgress(0, "Computing frequency.");
setProgress(25);
fq = frequency(model.param().evaluate("Lp"));
setProgress(75);
}
M. Hashiguchi and D. Mi 2018 124
セカント法（続き）
with(model.result("pg1"));
set("looplevel", new String[]{"7"}); // The first real eigenfrequency always the 7th computed eigenfrequency in solid mechanics
endwith();
model.result().numerical("gev1").setResult();
useGraphics(model.result(“pg1”), “graphics1”);
ウィンドウ(graphics1)
への表示zoomExtents(“graphics1”);
長さが変わったので画面にフィットさせるsetProgress(100);
プロセス100%
とするplay_sound();
音を出すmodel.setLastComputationTime(timeStamp()-t0); // record computation time solution_state = "solutionexists";
closeProgress();
プロットグループ１への反映
参考
COMSOL
のVideo
教材やBLOG
はぜひご覧ください。M. Hashiguchi and D. Mi 2018 126