把梦想照进现实

基于Matlab Optibox工具箱进行常微分方程参数拟合

2018.04.25

之前说了基于Berkeley的常微分方程参数拟合,由于该程序属于付费程序,可以进行计算,但是无法导出和保存。
本文介绍使用Matlab的Optibox工具箱实现常微分方程组系数拟合。

前提

Matlab2011 a或更高
Visual C++ 2015

安装

  • 安装Visual C++ 2015
  • 将Optibox解压缩于Matlab的工作目录,运行opyi_Install.m

帮助文件

http://www.i2c2.aut.ac.nz/Wiki/OPTI/index.php/Dynamic/DNLS

运行代码

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
% ODE System 常微分方程组
ode = @(t,z,p) [-p(1)*z(1) + 4;
2*z(1)  p(1)*z(2) + 5;
-4*z(1)  2*z(2)*z(3)  p(2)];
% Initial Conditions  初始值
z0 = [-1.5;1.25;1];
% True Parameter Values
p = [2.345;1.1];
% Generate Fitting Data
tm = 0:0.1:2; %measurement times  试验时间
odeInt = @(t,z) ode(t,z,p); %ODE function for ODE45
[~,zm] = ode45(odeInt,tm,z0); %Solve ODEs
% Starting Guess
theta0 = [1;0.5];  %%需要拟合的反应参数,初始值
% Create OPTI Object   最后拟合得到的反应参数为thetao
Opt = opti(ode,ode,data,tm,zm,z0,z0,theta0,theta0)
% Solve
[theta,fval,exitflag,info] = solve(Opt)
% Plot the Solution
plot(Opt)

我们需要用的

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
% ODE System 常微分方程组
ode = @(t,z,p) [-p(1)*z(1) + 4;
2*z(1)  p(1)*z(2) + 5;
-4*z(1)  2*z(2)*z(3)  p(2)];
% Initial Conditions 初始值
z0 = [-1.5;1.25;1];
% True Parameter Values
%p = [2.345;1.1];
% Generate Fitting Data
tm = 0:0.1:2; %measurement times 试验时间,本案例一共21个时间点
%odeInt = @(t,z) ode(t,z,p); %ODE function for ODE45
%[~,zm] = ode45(odeInt,tm,z0); %Solve ODEs
%zm为实验值
zm=[-1.50000000000000 1.25000000000000 1
-0.829847410728960 1.23122934889999 1.08891464236354
-0.299843823829989 1.32263477170971 0.937744166230725
0.119396419931630 1.47869265266965 0.638307388572418
0.450984569258225 1.66849705490454 0.268780940878075
0.713279518603320 1.87101913060957 -0.103393090731714
0.920731587587370 2.07274070551578 -0.433186879593826
1.08483205183455 2.26507683398946 -0.696928804994375
1.21462684427427 2.44317707294072 -0.890714278023006
1.31729157968141 2.60457626020426 -1.02217179283428
1.39849049014987 2.74849253183034 -1.10436440774115
1.46272077994400 2.87515877778552 -1.15117309475677
1.51352137402991 2.98551576019452 -1.17408059909133
1.55370586494275 3.08083451942504 -1.18228510531201
1.58548831552729 3.16258917236617 -1.18174229640595
1.61062899782998 3.23227896817056 -1.17669849596675
1.63051309925309 3.29138018651986 -1.16950533866256
1.64624190139260 3.34127127839269 -1.16170862782906
1.65868202108547 3.38322246677541 -1.15405079058414
1.66852245476047 3.41837172808451 -1.14696290618509
1.67630513208420 3.44773057475391 -1.14060057697597];
% Starting Guess
theta0 = [1;0.5]; %%需要拟合的反应参数,初始值
% Create OPTI Object 最后拟合得到的反应参数为thetao
Opt = opti(ode,ode,data,tm,zm,z0,z0,theta0,theta0)
% Solve
[theta,fval,exitflag,info] = solve(Opt)
% Plot the Solution
plot(Opt)

最终成图并得到结果,结果即为theta。

Comments
Write a Comment