MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

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--


  • Prev by Date: Re: Re: Delete problem
  • Next by Date: Re: Re: Integrate in version 5.1
  • Previous by thread: Re: Plotting a function of an interpolated function
  • Next by thread: Mathematica language issues