Struggling with FindFit for a model which takes results from coupled ODEs as arguments
- To: mathgroup at smc.vnet.net
- Subject: [mg108156] Struggling with FindFit for a model which takes results from coupled ODEs as arguments
- From: "Steven Mullen" <Steve at 21cm.com>
- Date: Tue, 9 Mar 2010 06:24:52 -0500 (EST)
I have been struggling for a few days to figure out how to estimate parameters in coupled ODEs when the results from the ODEs are used as input for the model in FindFit. It is a relatively straightforward mass transfer problem, with 2 components (NinWater and NinCPA1). These equations are designed to calculate the number of moles of substances per time, but I have data for total volume, so I am evaluating an equation to calculate the volume from (moles*molar volume). I have been trying to use FindFit (for a few days now) to estimate the parameters Lp and PCPA1 found in the ODEs to obtain a best fit to the volume data (with no luck). Here is the code I have so far. Lp := 2.64*^-27; (* mol^2/(min atm um^5) *) Rad = 63; (* um *) A = 4 p Rad^2; (* um^2*) R=0.08206; (* L atm /(mol K) *) T=286; minutes = 10; Vb = 0.19; Vbv = 4 /3 p Rad^3 * Vb (* um^3 *) volconv = 1*^15 (* converts um^3 to L ; 1x10^15 um^3 per L *); molconv = 1*^6 (* converts moles to micromoles ; 1x10^6 umol per mol *); PMVWater = 1.8*^7; (* um^3/umol*) PMVDMSO = 7.13*^7; PCPA1 := 4.38*^-31; BSalt = 2.749 (*salt*); BCPA1 = 2.423 (*DMSO*); insecondvirialcoefficientlist={BSalt, BCPA1} ; exsecondvirialcoefficientlist={BSalt, BCPA1} NinWater0 = 4/3 p Rad^3 (1-Vb)/PMVWater NinSalt = 1.663*2.062*^-4 (* Salts; assumes 0.3 osmolal =8A 0.176 molal *); NinCPA10 =0.1*^-100(* DMSO *); totmolin[t_] := (NinWater[t]+ NinSalt+ NinCPA1[t]) XinSalt[t_] := NinSalt/totmolin[t] XinCPA1[t_] := NinCPA1[t]/totmolin[t] insolutemolefxlist[t_] := {XinSalt[t],XinCPA1[t] waterin[t_]:= Total[insolutemolefxlist[t]] (1+insecondvirialcoefficientlist.insolutemolefxlist[t]) NexWater = 0.056*^6 (* water *); NexSalt = 1.663*245; (* assumes mole fraction of 0.003 *) NexCPA1 = 1694 (* DMSO *); totmolex = NexWater+ NexSalt+ NexCPA1; XexWater =NexWater/totmolex; XexSalt = NexSalt/totmolex; XexCPA1 = NexCPA1/totmolex; exsolutemolefxlist = {XexSalt,XexCPA1} waterex=Total[exsolutemolefxlist] * (1+exsecondvirialcoefficientlist.exsolutemolefxlist); V0=NinWater0 * PMVWater + Vbv + NinCPA10*PMVDMSO tref=T+273.16; volconv = 1*^15 (* ******************************************************************************************************************* *) (* The 2 parameters I wish to find are Lp and PCPA1 here in the 2 ODEs *) solution = NDSolve[{NinWater'[t]=8A- Lp*volconv*molconv*A*R*T (waterex - waterin[t]), NinCPA1'[t] =8A PCPA1*volconv*molconv*A*R*T (Log[XexCPA1]+ (0.5-BCPA1) XexWater (1-XexCPA1)- (0.5-BSalt) XexWater XexSalt- Log[(NinCPA1[t]/(NinWater[t]+NinSalt+NinCPA1[t]))]- (0.5-BCPA1) (NinWater[t]/(NinWater[t]+NinSalt+NinCPA1[t])) (1-(NinCPA1[t]/(NinWater[t]+NinSalt+NinCPA1[t])))+ (0.5-BSalt) (NinWater[t]/(NinWater[t]+NinSalt+NinCPA1[t])) (NinSalt/(NinWater[t]+NinSalt+NinCPA1[t]))), NinWater[0] =8ANinWater0, NinCPA1 [0] =8A NinCPA10},{NinWater, NinCPA1}, {t,0,minutes}] (* End of ODEs *) (*Start of Curve Fitting *) Vtotal[t] := (Vbv + NinWater[t]*PMVWater + NinCPA1[t]*PMVDMSO) /.solution; (* So I am making this calculation using output from the ODEs to determine the total volume, which includes a fixed volume (Vbv), the result from NinWater'[t] * molar volume of water, and likewise for the DMSO. The sum of all 3 components is the total volume, which is what I am trying to fit to the data (see below). volume:=Evaluate[Vtotal[t]] Clear[Lp,PCPA1] (* It seems that I have to clear these values here to make FindFit work to the extent it is working now) *) data = {1.04739=D7106,1.02856=D7106,1.01005=D7106,991863.,974012.,956510. ,939369.,922600.,906215.,890226.,874644.,859480.,844747.,830455.,816614.,803234.,790325.,777895.,765951.,754500.,743546.,733094., 723147.,713704.,704765.,696328.,688389.,680941.,673978.,667489.,661465.,655 892.,650759.,646048.,641746.,637835.,634299.,631118.,628276.,625754.,623534.,621597.,619927.,618505.,617316.,616342., 615569.,614981.,614565.,614307.,614195.,614216.,614360.,614616.,614975.,615427.,615965.,616581.,617268.,618019.,618827.,619689.,620598.,621549.,622540.,623564.,624620.,625704.,626812.,627942., 629092.,630259.,631441.,632636.,633842.,635059.,636284.,637516.,638754.,639997.,641244.,642494.,643746.,644999.,646253.,647508.,648762.,650015.,651267.,652517.,653765.,655011.,656254.,657494., 658730.,659963.,661193.,662419.,663641.,664858.,666072.,667281.,668485.,669686.,670881.,672072.,673258.,674440.,675617.,676789.,677956.,679118.,680275.,681428.,682575.,683718.,684856.,685989., 687117.,688240.,689359.,690472.,691581.,692685.,693784.,694879.,695969.,697054.,698134.,699210.,700281.,701347.,702409.,703467.,704519.,705568.,706612.,707651.,708686.,709717.,710743.,711765., 712783.,713797.,714806.,715811.,716812.,717809.,718801.,719790.,720774.,721755.,722731.,723704.,724673.,725637.,726598.,727555.,728508.,729458.,730403.,731345.,732283.,733218.,734148.,735076., 735999.,736919.,737835.,738748.,739657.,740563.,741466.,742365.,743260.,744152.,745041.,745926.,746808.,747687.,748562.,749435.,750304.,751169.,752032.,752891.,753748.,754601.,755451.,756298., 77142.,757982.,758820.,759655.,760487.,761315.,762141.,762964.,763784.,764601.,765415.,766227.,767035.,767841.,768644.,769444.,770241.,771036.,771828. ,772617.,773403.,774187.,774968.,775746., 776522.,777295.,778066.,778834.,779599.,780362.,781122.,781880.,782635.,783388.,784138.,784885.,785631.,786374.,787114.,787852.,788587.,789321.,790051.,790780.,791506.,792230.,792951.,793670., 794387.,795102.,795814.,796524.,797232.,797937.,798640.,799341.,800040.,800737.,801432.,802124.,802814.,803502.,804188.,804872.,805554.,806233.,806911.,807586.,808260.,808931.,809600.,810268., 810933.,811596.,812257.,812917.,813574.,814229.,814882.,815534.,816183.,816831.,817476.,818120.,818762.,819402.,820040.,820676.,821310.,821943.,822573.,823202.,823829.,824454.,825077.,825699., 826319.,826937.,827553.,828167.,828780.,829391.,830000.,830607.,831213.,831817.,832419.,833019.,833618.,834215.,834811.,835405.,835997.,836587.,837176.,837763.,838349.,838933.,839515.,840096., 840675.,841253.,841829.,842403.,842976.,843547.,844117.,844685.,845252.,845817.,846380.,846942.,847503.,848061.,848619.,849175.,849729.,850282.,850834.,851384.,851932.,852480.,853025.,853569., 854112.,854653.,855193.,855732.,856269.,856804.,857339.,857872.,858403.,858933.,859462.,859989.,860515.,861040.,861563.,862085.,862605.,863124.,863642.,864159.,864674.,865188.,865700.,866211., 866721.,867230.,867737.,868243.,868748.,869252.,869754.,870255.,870754.,871253.,871750.,872246.,872741.,873234.,873726.,874217.,874707.,875196.,875683.,876169.,876654.,877138.,877620.,878101., 878582.,879061.,879538.,880015.,880490.,880965.,881438.,881910.,882381.,882850.,883319.,883786.,884252.,884718.,885182.,885645.,886106.,886567.,887026.,887485.,887942.,888399.,888854.,889308., 889761.,890213.,890664.,891113.,891562.,892010.,892456.,892902.,893346.,893790.,894232.,894673.,895114.,895553.,895991.,896428.,896864.,897300.,897734.,898167.,898599.,899030.,899460.,899889., 900318.,900745.,901171.,901596.,902020.,902443.,902866.,903287.,903707.,904127.,904545.,904962.,905379.,905794.,906209.,906623.,907035.,907447.,907858.,908268.,908677.,909085.,909492.,909898., 910303.,910708.,911111.,911514.,911916.,912316.,912716.,913115.,913513.,913910.,914307.,914702.,915097.,915490.,915883.,916275.,916666.,917056.,917446.,917834.,918222.,918609.,918994.,919380., 919764.,920147.,920530.,920911.,921292.,921672.,922052.,922430.,922807.,923184.,923560.,923935.,924310.,924683.,925056.,925428.,925799.,926169.,926538.,926907.,927275.,927642.,928008.,928374., 928739.,929103.,929466.,929828.,930190.,930551.,930911.,931270.,931629.,931986.,932343.,932700.,933055.,933410.,933764.,934117.,934470.,934822.,935173.,935523.,935873.,936221.,936570.,936917., 937264.,937610.,937955.,938299.,938643.,938986.,939328.,939670.,940011.,940351.,940691.,941029.,941368.,941705.,942042.,942378.,942713.,943048.,943382.,943715.,944047.,944379.,944710.,945041., 945371.,945700.,946028.,946356.,946683.,947010.,947336.,947661.,947985.,948309.,948632.,948955.,949277.,949598.,949919.,950239.,950558.,950876.,951194.,951512.,951828.,952145.,952460.,952775., 953089.,953402.,953715.,954028.,954339.,954650.,954961.,955271.,955580.,955888.,956196.,956504.,956810.,957116.,957422.,957727.,958031.,958335.,958638.,958940.,959242.,959544.,959844.,960144., 960444.,960743.,961041.} (* These data represent one datapoint per second, for 10 minutes *) FindFit[data,volume,{{Lp,2.64*^-27},{PCPA1, 4.38*^-31}},t] (* When I evaluate FindFit, I get errors suggesting it is calculating beyond 10 minutes. *) I have been trying to use the online help and past posts, but with no luck. If anyone can assist me I would be very grateful. Thank you. Steven F. Mullen Ph.D.