lsqcurvefit to estimate parameters in ODE45

3 views (last 30 days)
Antonella Treglia
Antonella Treglia on 27 May 2021
Commented: Torsten on 27 May 2021
Hello!
I am trying to solve a parametric ODE system. I woul like to estimate the parameters by fitting the calculated data with experimental data.
Model:
dydt(1) = - theta(2).*y(1).*y(2)-theta(3).*y(1).*(theta(4)-y(3))-theta(5)*y(1).*y(2).^2;
dydt(2) = -theta(2).*y(1).*y(2)-theta(5)*y(1).*y(2).^2;
dydt(3) = theta(3).*y(1).*(theta(4)-y(3));
dydt(4) = 0;
it is a system of differential equations with 5 parameters (theta). These are constant parameters that affect the time evolution of the set of solutions y.
What I can empirically measure is the quantity
C=theta(2)*y(1)*y(2)/max(theta(2)*y(1)*y(2));
Therefore I would like to solve the system of ODE and fit C to the measured values to obtain the best guess for the set of parameters theta.
Can you please help me??
Here my code:
% for ODES calculation
tspan = [0 4e-6]; %seconds
precision = 0.00001;
%parameters initial values
theta0 = [4e17 2e-11 7e-8 5e15 9.3e-30];
% Variables
n_abs =3.3e15;
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
[Cfit,T] = kinetics(theta,tspan,precision);
function [C,T] = kinetics(theta,tspan,precision)
fprintf(1, '\t\tTheta(%d) = %8.5f\n',1, theta(1))
y0 = [3.3e15 3.3e15+theta(1) 0 0 ]; %initial conditions ne0, nh0, nte0, nth0
Opt = odeset('Events', @myEvent);
[T,y] = ode45(@(t,y) odefcn, tspan, y0, Opt);
function dydt = odefcn(t,y)
fprintf(1, '\t\tTheta(%d) = %8.5f\n',3, theta(3))
dydt = zeros(4,1);
dydt(1) = - theta(2).*y(1).*y(2)-theta(3).*y(1).*(theta(4)-y(3))-theta(5)*y(1).*y(2).^2;
dydt(2) = -theta(2).*y(1).*y(2)-theta(5)*y(1).*y(2).^2;
dydt(3) = theta(3).*y(1).*(theta(4)-y(3));
dydt(4) = 0;
end
function [value, isterminal, direction] = myEvent(t, y)
value = (y(1).*y(2) < precision*y0(1)*y0(2));
isterminal = 1; % Stop the integration
direction = 0;
end
C=theta(2)*y(1)*y(2)/max(theta(2)*y(1)*y(2));
end
I get this error
Not enough input arguments.
Error in x_solveODE>kinetics/odefcn (line 150)
dydt(1) = - theta(2).*y(1).*y(2)-theta(3).*y(1).*(theta(4)-y(3))-theta(5)*y(1).*y(2).^2;
Error in x_solveODE>@(t,y)odefcn (line 134)
[T,y] = ode45(@(t,y) odefcn, tspan, y0, Opt);
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in x_solveODE>kinetics (line 134)
[T,y] = ode45(@(t,y) odefcn, tspan, y0, Opt);
Error in lsqcurvefit (line 225)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in x_solveODE (line 113)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
Here the experimental data vectors:
t=[-0.341796825000000
-0.292968725000000
-0.244140625000000
-0.195312525000000
-0.146484325000000
-0.0976562249999997
-0.0488281250000000
-2.50000002921524e-08
0.0488281749999997
0.0976562750000003
0.146484375000000
0.195312475000000
0.244140675000000
0.292968775000000
0.341796875000000
0.390624975000000
0.439453175000000
0.488281275000000
0.537109375000000
0.585937475000000
0.634765675000000
0.683593775000000
0.732421875000000
0.781249975000000
0.830078175000000
0.878906275000000
0.927734375000000
0.976562475000000
1.02539067500000
1.07421877500000
1.12304687500000
1.17187497500000
1.22070317500000
1.26953127500000
1.31835937500000
1.36718747500000
1.41601567500000
1.46484377500000
1.51367187500000
1.56249997500000
1.61132817500000
1.66015627500000
1.70898437500000
1.75781247500000
1.80664067500000
1.85546877500000
1.90429687500000
1.95312497500000
2.00195317500000
2.05078127500000
2.09960937500000
2.14843747500000
2.19726567500000
2.24609377500000
2.29492187500000
2.34374997500000
2.39257817500000
2.44140627500000
2.49023437500000
2.53906247500000
2.58789067500000
2.63671877500000
2.68554687500000
2.73437497500000
2.78320317500000
2.83203127500000
2.88085937500000
2.92968747500000
2.97851567500000
3.02734377500000
3.07617187500000
3.12499997500000
3.17382817500000
3.22265627500000
3.27148437500000
3.32031247500000
3.36914067500000
3.41796877500000
3.46679687500000
3.51562497500000
3.56445317500000
3.61328127500000
3.66210937500000
3.71093747500000
3.75976567500000
3.80859377500000
3.85742187500000
3.90624997500000
3.95507817500000
4.00390627500000
4.05273437500000
4.10156247500000
4.15039067500000
4.19921877500000
4.24804687500000
4.29687497500000
4.34570317500000
4.39453127500000
4.44335937500000
4.49218747500000
4.54101567500000
4.58984377500000
4.63867187500000
4.68749997500000
4.73632817500000
4.78515627500000
4.83398437500000
4.88281287500000
4.93164087500000
4.98046887500000
5.02929687500000
5.07812487500000
5.12695287500000
5.17578087500000
5.22460987500000
5.27343787500000
5.32226587500000
5.37109387500000
5.41992187500000
5.46874687500000
5.51757787500000
5.56640587500000
5.61523487500000
5.66406287500000
5.71289087500000
5.76171887500000
5.81054687500000
5.85937487500000
5.90820287500000
5.95703087500000
6.00585987500000
6.05468787500000
6.10351587500000
6.15234687500000
6.20117187500000
6.24999987500000
6.29882787500000
6.34765587500000
6.39648487500000
6.44531287500000
6.49414087500000
6.54296887500000
6.59179687500000
6.64062487500000
6.68945287500000
6.73828087500000
6.78710987500000
6.83593787500000
6.88476587500000
6.93359387500000
6.98242187500000
7.03124687500000
7.08007787500000
7.12890587500000
7.17773487500000
7.22656287500000
7.27539087500000
7.32421887500000
7.37304687500000
7.42187487500000
7.47070287500000
7.51953087500000
7.56835987500000
7.61718787500000
7.66601587500000
7.71484687500000
7.76367187500000
7.81249987500000
7.86132787500000
7.91015587500000
7.95898487500000
8.00781287500000
8.05664087500000
8.10546887500000
8.15429687500000
8.20312487500000
8.25195287500000
8.30078087500000
8.34960987500000
8.39843787500000
8.44726587500000
8.49609387500000
8.54492187500000
8.59374687500000
8.64257787500000
8.69140587500000
8.74023487500000
8.78906287500000
8.83789087500000
8.88671887500000
8.93554687500000
8.98437487500000
9.03320287500000
9.08203087500000
9.13085987500000
9.17968787500000
9.22851587500000
9.27734687500000
9.32617187500000
9.37499987500000
9.42382787500000
9.47265587500000
9.52148487500000
9.57031287500000
9.61914087500000
9.66796887500000
9.71679687500000
9.76562487500000
9.81445287500000
9.86328087500000
9.91210987500000
9.96093787500000
10.0097658750000
10.0585938750000
10.1074218750000
10.1562468750000
10.2050778750000
10.2539058750000
10.3027348750000
10.3515628750000
10.4003908750000
10.4492188750000
10.4980468750000
10.5468748750000
10.5957028750000
10.6445308750000
10.6933598750000
10.7421878750000
10.7910158750000
10.8398468750000
10.8886718750000
10.9374998750000
10.9863278750000
11.0351558750000
11.0839848750000
11.1328128750000
11.1816408750000
11.2304688750000
11.2792968750000
11.3281248750000
11.3769528750000
11.4257808750000
11.4746098750000
11.5234378750000
11.5722658750000
11.6210938750000
11.6699218750000
11.7187468750000
11.7675778750000
11.8164058750000
11.8652348750000
11.9140628750000
11.9628908750000
12.0117188750000
12.0605468750000
12.1093748750000
12.1582028750000
12.2070308750000
12.2558598750000
12.3046878750000
12.3535158750000
12.4023468750000
12.4511718750000
12.4999998750000
12.5488278750000
12.5976558750000
12.6464848750000
12.6953128750000
12.7441408750000
12.7929688750000
12.8417968750000
12.8906248750000
12.9394528750000
12.9882808750000
13.0371098750000
13.0859378750000
13.1347658750000
13.1835938750000
13.2324218750000
13.2812468750000
13.3300778750000
13.3789058750000
13.4277348750000
13.4765628750000
13.5253908750000
13.5742188750000
13.6230468750000
13.6718748750000
13.7207028750000
13.7695308750000
13.8183598750000
13.8671878750000
13.9160158750000
13.9648468750000
14.0136718750000
14.0624998750000
14.1113278750000
14.1601558750000
14.2089848750000
14.2578128750000
14.3066408750000
14.3554688750000
14.4042968750000
14.4531248750000
14.5019528750000
14.5507808750000
14.5996098750000
14.6484378750000
14.6972658750000
14.7460938750000
14.7949218750000
14.8437468750000
14.8925778750000
14.9414058750000
14.9902348750000
15.0390628750000
15.0878908750000
15.1367188750000
15.1855468750000
15.2343748750000
15.2832028750000
15.3320308750000
15.3808598750000
15.4296878750000
15.4785158750000
15.5273468750000
15.5761718750000
15.6249998750000
15.6738278750000
15.7226558750000
15.7714848750000
15.8203128750000
15.8691408750000
15.9179688750000
15.9667968750000
16.0156248750000
16.0644528750000
16.1132808750000
16.1621098750000
16.2109378750000
16.2597658750000
16.3085938750000
16.3574218750000
16.4062468750000
16.4550778750000
16.5039058750000
16.5527348750000
16.6015628750000
16.6503908750000
16.6992188750000
16.7480468750000
16.7968748750000
16.8457028750000
16.8945308750000
16.9433598750000
16.9921878750000
17.0410158750000
17.0898468750000
17.1386718750000
17.1874998750000
17.2363278750000
17.2851558750000
17.3339848750000
17.3828128750000
17.4316408750000
17.4804688750000
17.5292968750000
17.5781248750000
17.6269528750000
17.6757808750000
17.7246098750000
17.7734378750000
17.8222658750000
17.8710938750000
17.9199218750000
17.9687468750000
18.0175778750000
18.0664058750000
18.1152348750000
18.1640628750000
18.2128908750000
18.2617188750000
18.3105468750000
18.3593748750000
18.4082028750000
18.4570308750000
18.5058598750000
18.5546878750000
18.6035158750000
18.6523468750000
18.7011718750000
18.7499998750000
18.7988278750000
18.8476558750000
18.8964848750000
18.9453128750000
18.9941408750000
19.0429688750000
19.0917968750000
19.1406248750000
19.1894528750000
19.2382808750000
19.2871098750000
19.3359378750000
19.3847658750000
19.4335938750000
19.4824218750000
19.5312468750000
19.5800778750000
19.6289058750000
19.6777348750000
19.7265628750000
19.7753908750000
19.8242188750000
19.8730468750000
19.9218748750000
19.9707028750000
20.0195308750000
20.0683598750000
20.1171878750000
20.1660158750000
20.2148468750000
20.2636718750000
20.3124998750000
20.3613278750000
20.4101558750000
20.4589848750000
20.5078128750000
20.5566408750000
20.6054688750000
20.6542968750000
20.7031248750000
20.7519528750000
20.8007808750000
20.8496098750000
20.8984378750000
20.9472658750000
20.9960938750000
21.0449218750000
21.0937468750000
21.1425778750000
21.1914058750000
21.2402348750000
21.2890628750000
21.3378908750000
21.3867188750000
21.4355468750000
21.4843748750000
21.5332028750000
21.5820308750000
21.6308598750000
21.6796878750000
21.7285158750000
21.7773468750000
21.8261718750000
21.8749998750000
21.9238278750000
21.9726558750000
22.0214848750000
22.0703128750000
22.1191408750000
22.1679688750000
22.2167968750000
22.2656248750000
22.3144528750000
22.3632808750000
22.4121098750000
22.4609378750000
22.5097658750000
22.5585938750000
22.6074218750000
22.6562468750000]; %experimental data (time axis)
c=[0.918600438554676
0.962551244160549
0.951301363333016
0.996110210696921
0.974945180665459
0.988864524740204
0.982000190675946
0.997254266374297
0.978949375536276
0.978377347697588
0.964744017542187
0.943006959672037
0.941386214129088
0.932043092763848
0.927943559919916
0.900390885689770
0.910782724759272
0.904109066641243
0.893240537706168
0.889427018781581
0.877986462007818
0.838421203165221
0.860348936981600
0.847573648584231
0.804576222709505
0.810391839069501
0.794279721613119
0.786462007817714
0.786366669844599
0.780455715511488
0.777881590237392
0.748898846410525
0.726971112594146
0.722108875965297
0.706854800266946
0.702278577557441
0.679874153875489
0.654514253026981
0.666717513585661
0.644980455715511
0.634397940699781
0.628963676232243
0.620287920678806
0.589589093335876
0.596358089427019
0.584345504814568
0.566040613976547
0.557936886261798
0.555934788826390
0.549261130708361
0.556220802745734
0.548307750977214
0.525998665268376
0.532195633520831
0.523043188101821
0.508647154161503
0.512842024978549
0.516655543903137
0.514176756602155
0.513700066736581
0.493297740490037
0.499018018876919
0.488721517780532
0.478234340737916
0.475469539517590
0.470321288969397
0.473467442082181
0.454304509486128
0.463266278958909
0.450300314615311
0.446010105825150
0.438764419868434
0.429897988368767
0.438573743922204
0.421794260654019
0.419982839164839
0.420745542949757
0.424177709981886
0.403870721708456
0.399675850891410
0.392334826961579
0.392906854800267
0.383945085327486
0.375555343693393
0.370407093145200
0.362398703403566
0.362112689484222
0.354676327581276
0.353055582038326
0.326646963485556
0.345142530269806
0.325502907808180
0.326837639431786
0.317685194012775
0.307102678997045
0.299952331013443
0.301763752502622
0.301668414529507
0.295757460196396
0.281742778148537
0.281075412336734
0.274306416245591
0.271732290971494
0.263533225283630
0.268490799885594
0.258289636762322
0.265821336638383
0.274401754218705
0.255715511488226
0.256192201353799
0.252378682429212
0.247516445800362
0.246086376203642
0.243798264848889
0.242082181332825
0.249613881208885
0.238364000381352
0.236457240919058
0.235503861187911
0.234359805510535
0.229783582801030
0.228067499284965
0.228639527123653
0.227495471446277
0.224444656306607
0.226446753742015
0.216817618457432
0.213290113452188
0.212336733721041
0.217866336161693
0.201849556678425
0.202326246543998
0.207474497092192
0.196510630184002
0.209095242635142
0.191457717608924
0.195938602345314
0.190027648012203
0.186690818953189
0.183926017732863
0.179254457050243
0.180493850700734
0.167909238249595
0.169911335685003
0.169339307846315
0.162951663647631
0.164763085136810
0.163523691486319
0.162951663647631
0.158375440938126
0.149699685384689
0.147983601868624
0.143407379159119
0.142549337401087
0.144074744970922
0.139784536180761
0.143121365239775
0.139307846315187
0.129869386976833
0.140547239965678
0.132824864143388
0.134064257793879
0.133206216035847
0.125865192106016
0.126341881971589
0.124625798455525
0.120907617504052
0.123100390885690
0.122623701020116
0.123100390885690
0.121956335208313
0.123481742778149
0.120430927638478
0.116903422633235
0.117856802364382
0.117952140337496
0.116808084660120
0.114424635332253
0.116903422633235
0.114519973305368
0.112994565735532
0.115854704928973
0.111469158165697
0.108609018972257
0.104986175993898
0.105653541805701
0.103842120316522
0.100028601391934
0.101840022881114
0.0988845457145581
0.0917341977309562
0.103746782343407
0.0907808179998093
0.0867766231289923
0.0883020306988274
0.0895414243493183
0.0883973686719420
0.0863952712365335
0.0884927066450567
0.0822004004194871
0.0840118219086662
0.0839164839355515
0.0753360663552293
0.0754314043283440
0.0775288397368672
0.0763847840594909
0.0732386309467061
0.0709505195919535
0.0719992372962151
0.0697111259414625
0.0679950424253980
0.0692344360758890
0.0689484221565450
0.0663742968824483
0.0683763943178568
0.0704738297263800
0.0677090285060540
0.0653255791781867
0.0663742968824483
0.0674230145867099
0.0651349032319573
0.0633234817427782
0.0661836209362189
0.0659929449899895
0.0630374678234341
0.0646582133663838
0.0612260463342549
0.0630374678234341
0.0648488893126132
0.0594146248450758
0.0633234817427782
0.0615120602535990
0.0552197540280294
0.0592239488988464
0.0611307083611403
0.0577938793021260
0.0582705691676995
0.0551244160549147
0.0533129945657355
0.0494994756411479
0.0578892172752407
0.0551244160549147
0.0503575173991801
0.0508342072647536
0.0525502907808180
0.0579845552483554
0.0518829249690152
0.0510248832109829
0.0506435313185242
0.0476880541519687
0.0534083325388502
0.0478787300981981
0.0519782629421299
0.0477833921250834
0.0484507579368863
0.0480694060444275
0.0504528553722948
0.0490227857755744
0.0503575173991801
0.0460673086090190
0.0478787300981981
0.0472113642863953
0.0448279149585280
0.0480694060444275
0.0453999427972161
0.0491181237486891
0.0416817618457432
0.0429211554962342
0.0441605491467251
0.0435885213080370
0.0438745352273811
0.0439698732004958
0.0440652111736104
0.0435885213080370
0.0435885213080370
0.0445419010391839
0.0444465630660692
0.0409190580608256
0.0395843264372199
0.0361521594050911
0.0385356087329583
0.0410143960339403
0.0396796644103346
0.0397750023834493
0.0380589188673849
0.0372008771093527
0.0400610163027934
0.0396796644103346
0.0374868910286967
0.0326246543998475
0.0371055391362380
0.0358661454857470
0.0333873581847650
0.0351987796739441
0.0331966822385356
0.0354847935932882
0.0267137000667366
0.0323386404805034
0.0311945848031271
0.0310039088568977
0.0300505291257508
0.0289064734483745
0.0301458670988655
0.0288111354752598
0.0274764038516541
0.0290018114214892
0.0279530937172276
0.0254743064162456
0.0294785012870626
0.0229955191152636
0.0281437696634570
0.0261416722280484
0.0259509962818190
0.0255696443893603
0.0253789684431309
0.0237582229001811
0.0269997139860807
0.0261416722280484
0.0232815330346077
0.0212794355991992
0.0243302507388693
0.0220421393841167
0.0242349127657546
0.0219468014110020
0.0243302507388693
0.0245209266850987
0.0203260558680522
0.0190866622175613
0.0215654495185432
0.0222328153303461
0.0195633520831347
0.0195633520831347
0.0207074077605110
0.0198493660024788
0.0198493660024788
0.0217561254647726
0.0205167318142816
0.0185146343788731
0.0192773381637906
0.0200400419487082
0.0218514634378873
0.0179426065401850
0.0200400419487082
0.0192773381637906
0.0190866622175613
0.0133663838306798
0.0170845647821527
0.0154638192392030
0.0167985508628087
0.0171799027552674
0.0150824673467442
0.0134617218037945
0.0181332824864143
0.0153684812660883
0.0151778053198589
0.0169892268090380
0.0138430736962532
0.0144151015349414
0.0156544951854324
0.0154638192392030
0.0136523977500238
0.0147011154542854
0.0139384116693679
0.0121269901801888
0.0162265230241205
0.0120316522070741
0.0134617218037945
0.0147011154542854
0.0149871293736295
0.0113642863952712
0.0105062446372390
0.0149871293736295
0.0106015826103537
0.0134617218037945
0.0115549623415006
0.0131757078844504
0.0115549623415006
0.0124130040995328
0.0134617218037945
0.0115549623415006
0.0103155686910096
0.0112689484221565
0.0101248927447802
0.0118409762608447
0.0121269901801888
0.0119363142339594
0.0108875965296978
0.0102202307178949
0.0127943559919916
0.0121269901801888
0.00993421679855086
0.0118409762608447
0.00974354085232148
0.00888549909428926
0.00964820287920679
0.00983887882543617
0.00945752693297740
0.00898083706740395
0.00898083706740395
0.0105062446372390
0.00831347125560111
0.00831347125560111
0.00869482314805987
0.00669272571265135
0.00859948517494518
0.00821813328248642
0.00678806368576604
0.00802745733625703
0.00802745733625703
0.00669272571265135
0.00688340165888073
0.00488130422347221
0.00631137382019258
0.00593002192773382
0.00764610544379827
0.00707407760511012
0.00621603584707789
0.00669272571265135
0.00688340165888073
0.00669272571265135
0.00573934598150443
0.00526265611593098
0.00659738773953666
0.00640671179330728
0.00802745733625703
0.00583468395461913
0.00650204976642197
0.00573934598150443
0.00545333206216036
0.00640671179330728
0.00659738773953666
0.00678806368576604
0.00583468395461913
0.00688340165888073
0.00774144341691296
0.00612069787396320
0.00612069787396320
0.00659738773953666
0.00659738773953666
0.00621603584707789
0.00688340165888073
0.00669272571265135
0.00640671179330728
0.00707407760511012
0.00640671179330728
0.00478596625035752
0.00640671179330728
0.00793211936314234
0.00669272571265135
0.00526265611593098
0.00640671179330728
0.00621603584707789
0.00373724854609591
0.00383258651921060
0.00736009152445419
0.00402326246543999
0.00402326246543999
0.00364191057298122
0.00449995233101344
0.00459529030412813
0.00383258651921060
0.00373724854609591
0.00230717894937554
0.00488130422347221
0.00411860043855468
0.00488130422347221
0.00411860043855468
0.00402326246543999
0.00268853084183430
0.00478596625035752
0.00278386881494899
0.00326055868052245]; %experimental data (yaxis)

Accepted Answer

Jan
Jan on 27 May 2021
Modify the line
[T,y] = ode45(@(t,y) odefcn, tspan, y0, Opt);
to
[T,y] = ode45(@(t,y) odefcn(t, y), tspan, y0, Opt);
or the faster version:
[T,y] = ode45(@odefcn, tspan, y0, Opt);
With the original code, odefcn is called without input arguments.
  2 Comments
Jan
Jan on 27 May 2021
The function to be optimized is called with the 2 inputs "myfun(x,xdata)". You want to provide a 3rd input, so you have to define it:
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = ...
lsqcurvefit(@(x, xdata) kinetics(x, xdata, precision), theta0, t, c);
Matlab cannot guess, how many inputs an anonymous function needs. You have to declare this explicitly.

Sign in to comment.

More Answers (1)

Torsten
Torsten on 27 May 2021
Modifications in the following lines are necessary:
[ ... ] = lsqcurvefit(@(theta,t) kinetics(theta,t,precision),...)
function C = kinetics(theta,t,precision)
[T,Y] = ode15s(@odefun,t,y0,Opt)
C = theta(2)*Y(:,1).*Y(:,2)/max(theta(2)*Y(:,1).*Y(:,2))
  6 Comments
Torsten
Torsten on 27 May 2021
@Star Strider - Thank you for the tip - really a usful feature.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!