DyMat は Dymola や OpenModelica, JModelica.org などの結果ファイル(matファイル)を読み込んだり、他のファイル形式にエクスポートすることができるライブラリです。これを使って JModelica.org によるシミュレーション結果ファイルのデータをプロットします。DyMat の詳細な使用方法は、Dymat User Manual Dy Mat Guideに記載されています。
まず、mat ファイルのリストを表示します。
%ls *.mat
CSTR_CSTR_Init_Optimization_result.mat CSTR_CSTR_result.mat ClassExample1_MassA_result.mat ClassExample2_MassB_result.mat EngineV6_analytic_with_input_result.mat Modelica_Mechanics_Rotational_Examples_First_result.mat RLC_Circuit_result.mat VDP_result.mat
Van der Pole oscillator のシミュレーション結果 VDP_resullt.mat を読み込みます。
import DyMat
d = DyMat.DyMatFile('VDP_result.mat')
すべての変数名を表示します。
for name in d.names():
print(name)
_use_Brent_in_1d _rescale_after_singular_jac x2_0 x1_0 _nle_solver_use_last_integrator_step _block_solver_profiling _nle_active_bounds_mode _rescale_each_step _events_tol_factor _block_solver_experimental_mode _nle_solver_default_tol _nle_solver_check_jac_cond _use_newton_for_brent _cs_solver _block_jacobian_check_tol _cs_step_size _events_default_tol _block_jacobian_check _iteration_variable_scaling _nle_jacobian_update_mode _nle_solver_max_iter der(x1) _nle_solver_step_limit_factor _nle_solver_use_nominals_as_fallback _nle_solver_min_residual_scaling_factor _nle_jacobian_calculation_mode _nle_solver_exit_criterion _time_events_default_tol _residual_equation_scaling _nle_jacobian_finite_difference_delta _log_level x2 _nle_brent_ignore_error x1 _nle_solver_regularization_tolerance _nle_solver_max_residual_scaling_factor _nle_solver_min_tol _use_jacobian_equilibration _cs_experimental_mode _runtime_log_to_file u _enforce_bounds _nle_solver_tol_factor _cs_rel_tol der(x2) _nle_solver_max_iter_no_jacobian
アンダーバーから始まる変数名を除いて表示します。
for name in d.names():
if name.find('_') != 0:
print(name)
x2_0 x1_0 der(x1) x2 x1 u der(x2)
プロットするデータを抽出します。
t = d.abscissa('x1', valuesOnly=True)
x1 = d.data('x1')
x2 = d.data('x2')
データをプロットします。
%matplotlib notebook
import matplotlib.pyplot as plt
plt.figure(1)
plt.plot(t, x1, t, x2)
plt.legend(( 'x1', 'x2'))
plt.title('Van der Pol scillator.')
plt.ylabel('Angle (rad)')
plt.xlabel('Time (s)')
Text(0.5,0,u'Time (s)')