NonlinearRegress+NIntegrate
- To: mathgroup at smc.vnet.net
- Subject: [mg52945] NonlinearRegress+NIntegrate
- From: Ole Hölck <ole.hoelck at bam.de>
- Date: Fri, 17 Dec 2004 05:19:37 -0500 (EST)
- Reply-to: ole.hoelck at web.de
- Sender: owner-wri-mathgroup at wolfram.com
Hello,
i am having problems getting mathematica to fit a model function to
experimental data. the model is
function[p_,Vh0_]:= Nintegrate[gamma*(Vg-\
(Vg+3*(G-G0my0)/(4.0*mushear)+(Vg-Vh0)^2/(2.0*Vh0)\
-Sqrt[((Vg+3*(G-G0my0)/(4*mushear)+(Vg-Vh0)^2/(2*Vh0))^2)-Vg^2]))\
*1/(sigma*Sqrt[Pi])*\
(Exp[-((G-G0my0)/sigma)^2])/\
(1+Exp[(G-R*T*Log[p])/(R*T)]),\
{G,Glow,Ghi}]
Vh0 beeing the fit-parameter and p the independent variable.
i have already read of solutions concerning the use of NonlinearRegress
on a Nintegrate like using the option FindMinimum or teaching
mathematica the derivatives, but it seems they do not work here.
the errormessage is, that the integral is not numerical. however, if i
let mathematica print the function with the starting-value Vh0start i
get a nice curve through the data.
how do i trick mathematica to evaluate the integral while fitting?
i would be grateful if someone could help me.
best regards,
o. hoelck
--------------090708080009070000070401
filename="nonlinfit_nintegrate.matscript"
(*****************************************************************************
* *
* Fit of the Parameter Vh0 *
* *
* *
******************************************************************************)
Print["begin packages"](***********************************************************)
Remove["Global`*"]
Needs["Statistics`NonlinearFit`"]
Print["end packages - begin parameters"](******************************************)
NA=6.022*10^23 (* 1/mol *)
N0=6.7*10^21 (* 1/cm^3(Poly) *) (* # of holes/VolumePol. *)
Vmol=22.4*10^3 (* cm^3(Gas)/mol *) (* molar volume @ T=0 C, p=1 atm *)
R=8.31451 (* N m^3/mol K, Jm^2/mol K *)
p0=1.013 (* bar *) (* =1atm), hier ist my=my0 ! *)
mu0=0.0 (* chemical potential at reference-state, s. mamo25_kirchheim (1992)*)
gamma=1.5 (* =3(1-ny)/(1+ny), ny=poisson ratio=1/3 (glassy polymers) *)
Vg=46.0 (* cm^3/mol = gasmolecule volume (part.mol.vol. in liquid*)
Vh0start=15 (* cm^3/mol = starting value for Vh0-fit *)
mushear=1000/1.5 (* MPa, Shearmodulus of the Polymer *)
sigma=11600.0 (* J/mol, width of site-energy-distribution,result from fit_c_p_kir.m *)
G0my0=20000.0 (* J/mol, mean site energy, result from fit_c_p_kir.m *)
Glow=G0my0-2*mushear/3*(Vg-Vh0start)^2/(Vh0start)
Ghi=G0my0+2*mushear/3*(Vg-Vh0start)^2/(Vh0start)
T =273.15+35 (* K *)
pmin =1 (* bar *)
pmax =40 (* bar *)
Glow=G0my0-2*mushear/3*(Vg-Vh0start)^2/(Vh0start)
Ghi=G0my0+2*mushear/3*(Vg-Vh0start)^2/(Vh0start)
datafile="M:/mathematikOmat/dilatfit/experimental.dat"
resultfile="M:/mathematikOmat/dilatfit/result.dat"
Print["end parameters - begin function"](*****************************************)
(* deltaVgraf is used to plot the function with the starting value*)
deltaVgraf[pgraf_,Vh0graf_]:=NIntegrate[gamma*(Vg-\
(Vg+3*(G-G0my0)/(4.0*mushear)+(Vg-Vh0graf)^2/(2.0*Vh0graf)\
-Sqrt[((Vg+3*(G-G0my0)/(4*mushear)+(Vg-Vh0graf)^2/(2*Vh0graf))^2)-Vg^2]))\
*1/(sigma*Sqrt[Pi])*\
(Exp[-((G-G0my0)/sigma)^2])/\
(1+Exp[(G-R*T*Log[pgraf])/(R*T)]),\
{G,Glow,Ghi}]
(* this is the model function to be fitted to the datapoints *)
deltaV[p_,Vh0_]:=NIntegrate[gamma*(Vg-\
(Vg+3*(G-G0my0)/(4.0*mushear)+(Vg-Vh0)^2/(2.0*Vh0)\
-Sqrt[((Vg+3*(G-G0my0)/(4*mushear)+(Vg-Vh0)^2/(2*Vh0))^2)-Vg^2]))\
*1/(sigma*Sqrt[Pi])*\
(Exp[-((G-G0my0)/sigma)^2])/\
(1+Exp[(G-R*T*Log[pgraf])/(R*T)]),\
{G,Glow,Ghi}]
Print["end function - begin plot"](***********************************************)
(* read data *)
data=ReadList[datafile,{Number,Number}]
(* plot curve with starting value and data-points *)
startcurve=Plot[deltaVgraf[grafp,Vh0start], {grafp,pmin,pmax},
PlotStyle->Automatic,
DisplayFunction->Identity ]
points=ListPlot[data,DisplayFunction->Identity,
PlotStyle->{{PointSize[0.2]}} ]
Show[points,startcurve,DisplayFunction->$DisplayFunction]
Print["end plot - begin fit"](*****************************************************)
lsg=BestFitParameters/.NonlinearRegress[data,\
deltaV[p_,Vh0_],p,{Vh0,Vh0start},\
ShowProgress->True]
Print["end fit - begin plot"](*****************************************************)
(* plot resulting curve and datapoints *)
resultcurve=Plot[deltaV[p,Vh0]/.lsg, {p,pmin,pmax},
PlotStyle->Automatic,
DisplayFunction->Identity ]
Show[points,resultcurve,DisplayFunction->$DisplayFunction ]
pfile=1;
fh=OpenWrite[resultfile]
WriteString[fh,"fit ",datafile]; Write[fh]
Write[fh,lsg];Write[fh]
WriteString[fh," p C"]; Write[fh]
WriteString[fh," 0.0 0.0"]; Write[fh]
While[pfile<51,WriteString[fh,pfile," ",Evaluate[deltaV[pfile,Vh0]/.lsg]];Write[fh];++pfile]
Close[fh]
--------------090708080009070000070401
filename="experimental.dat"
MC44OAkwLjYxOTU4DQoxLjcJMC44Mjg3OA0KMi43CTEuMDYwOTMNCjMuNzEJMS4yNzEwNQ0K
NC41NQkxLjQ2MTc3DQo1LjkJMS42NjU3Ng0KNy45OQkxLjk5Mjc5DQo5Ljk1CTIuMjkyMTkN
CjE0LjM1CTIuOTAzNjENCjE3LjYJMy40MTg0OA0KMjAuNTkJMy42NzI4Nw0K
--------------090708080009070000070401--