常微分方程组系数拟合

2017-05-25 2017

之前采用Berkeley Madonna程序进行了常微分方程组系数拟合，由于该程序属于付费程序，可以进行计算，但是无法导出和保存，可以进行计算，并采用得到的反应速率K值，进行常微分方程组的计算。
 本文介绍使用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 

% 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)
  • 我们需要用的是 
 % 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)
Comments
Write a Comment
△