How do I solve two ode (dct/dt and dR/dt) simultaneously with ODE45 and optimize two parameters using global optimization ?

[darova]

Try ode45 and show your attempts

[Gokul Rao]

function dYdt = concentration(t,Y,p)

R=Y(1);

CT=Y(2);

%Known Parameters

r= 250; %0.684 um radius of adsorbent

rf= 0.684; %Radius of concentration front of solute penetrating adsorbent

P=2.439; %2.439 density of adsorbent

m=1; %1 g mass of adsorbent

g=50; %50 mL/min gas flowrate

CO=0.1; %0.1 initial CO2 concentration

qet=138.72; %138.72 CO2 concentration at equilibrium

kF=43.764; %43.764 Freundlich constant

nF=2.7697;

cet=24.416;

%Unknown Parameters

kg = p(1);

De = p(2);

%nonlinear differential equation to show change in concentration vs time

Ch= (qet*m)/(CO*g);

Bi= (kg/De)*r;

CET= cet/CO;

R= rf/r;

N= (3*Ch*qet*r^2)+((kF/nF)*(CET^(1/nF-1))*(Ch*Bi*(1-r^3)*CT)*(1/(r+Bi*(1-r))^2));

M= 1+Ch*(1-r^3)*(kF/nF)*(CET^(1/nF-1))*(Bi*(1-r))*(1/(r+Bi*(1-r)));

dRdt= -(Bi*(CO/P*qet)*(CT-CET))*(1/R^2);

dCTdt= N/M*(dRdt);

dYdt= [dCTdt;dRdt];

end

% Define the objective

function obj = objective(p)

%Initial condition at t=0

R=0.999;

CT=1.0;

Y0= [R;CT];

%Specific time points corresponding to experimental time points (minutes)

tspan= [0 5 10 15 20 25 30 35 40 45 50 55 60];

%function handle to pass p through substrate function and obtain numerical

%model by solving nonlinear differential equation

[t,Y]= ode45(@(t,Y)concentration(t,Y,p),tspan,Y0);

%Experimental data at each time points ts

CT1measured= [1 2 3 4 5 6 7 8 9 10 11 12 13];

CT1measured= CT1measured';

%Objective function to minimise the square of the difference between the

%numerical model and experimental data

A = (CT-CT1measured).^2;

obj = sum(A);

end

clear;

clc;

close all;

%Parameter Initial Guess

kg=1;

De=1;

p0= [kg;De];

%Objective function is the square of the difference between numerical and

fun = @objective; % function handle created

%identify unknown parameters

%[p,fval], find both the location of the minimum p of the nonlinear

%function and the value of the function at that minimum.

[p,fval] = fminunc(fun,p0);

%Optimised parameter values

kg = p(1);

De = p(2);

disp(['kg: ' num2str(kg)])

disp(['De: ' num2str(De)])

%Calculate model with updated parameters

p = [kg;De];

tspan= [0 5 10 15 20 25 30 35 40 45 50 55 60];

R=0.999;

CT=1.0;

Y0= [R;CT];

[t,Y]= ode45(@(t,Y)concentration(t,Y,p),tspan,Y0);

CT1measured= [1 2 3 4 5 6 7 8 9 10 11 12 13];

%Plot of updated parameter values

plot(tspan,CT1measured,'o',tspan,Y,'r');xlabel('time (min)');ylabel('Ct/Co');

title('Ct/Co vs time (min)');grid on

Answer 1

[Sulaymon Eshkabilov]

These are coupled ODEs. To solve these ODEs numerically:

(0) Specify all constant parameters

(1) Build fcn file or anonymous fcn handle

(2) Call and solve the fcn using one of these solvers: ode23, ode45, ode113

Answer 2

[Alan Weiss]

For a relevant example, see Optimize an ODE in Parallel.

Alan Weiss

MATLAB mathematical toolbox documentation

Answer 3

[Star Strider]

See for one example: Parameter Estimation for a System of Differential Equations and How to properly model a kinetic reaction system from experimental data using lsqcurvefit? for another.