nonlinearpiecewise.mws

---

An Algorithm For Fitting Nonlinear Piecewise Functions

Stanislav Barton              barton@mendelu.cz

Dept. of Physics             ++ 420 5 520844 - privat & fax

Mendel Univ.                 ++ 420 5 4513 2127 - office

Zemedelska 1

613 00 Brno  Czech Rep.                                                                 

This worksheet shows how to find parameters of two nonlinear functions fitting the given points . The first function   is valid on the range x  = <0.. >, the second one    on the range x  = < ). Where a1,a2,b1,b2 are linear and  c1, c2  are nonlinear parameters. Both functions have to satisfy continuity conditions at the point , ie.  and . These conditions enable us to determine the two parameters b1  and b2  as a functions  and . The remaining linear parameters a1, a2, c1, c2 and point of the continuity  could be computed using the Least Squares Method minimizing the squared sum of the residuals S , ie.  by solving its normal equations ie.: , , ,  and . The point  splits data vectors X and Y into two subvectors X1 and corresponding Y1 for , X2 and corresponding Y2 for . The equations Ea1  and Ea2  are linear in a1  and a2 , so we can solve them:  and . The remaining equations Ec1, Ec2  and Exi  are nonlinear for c1 , c2  and , so we shall solve them using Newton's method. To save space we shall use the following notation:  for the equation Sc1  and very simillar for Sc2 and Sxi . Newton's method is based on the linerization of functions, se we have to build the set of 3 equations:  

 

.

If we assign proper initial numerical values ,   and , we can solve these equations for   the corrections  dc1 , dc2  and dxii  and to build new initial values ,  and  and we can compute again new corrections. The procedure terminates when these corrections become smaller than the desired precision.

>	restart; with(linalg): with(plots): Maple start

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

>	g:=45+650*exp(-exp(-t/250)*6);  Widely used Gompertz function - Pattern of the growth of a bull. We shall try to approximate it.

>	plot(g,t=0..1500,title=`Gompertz function`,labels=["t [days]","M [kg]"]);

>	f1:=t->a1+b1*exp(t/c1); f2:=t->a2-b2*exp(-t/c2); Definition of new functions,   f1(t)  for    and f2(t)  for

>	e1:=f1(xi)=f2(xi);   e2:=diff(f1(xi),xi)=diff(f2(xi),xi); Continuity condition and continuity of the 1st derivatives

>	sol:=solve({e1,e2},{b1,b2}); Determination of two linear coefficients

>	assign(sol);

>	s1:=expand(K*(M1[i]-f1(T1[i]))^2); s2:=expand(K*(M2[i]-f2(T2[i]))^2);  Squared residuals of the 1st and 2nd functions. K  is necessary to create products at each element of the result

>	s1:=map(u->remove(has,u,i)           *simplify(sum(select(has,u,i),i=1..N1)),s1); s2:=map(u->remove(has,u,i)           *simplify(sum(select(has,u,i),i=1..N2)),s2);  Residuals are summed.

>	K:=1:  K is no longer necessary

>

DE:=[diff(E1(c1,xi),c1,c1)=E1c1c1, diff(E1(c1,xi),xi,xi)=E1xixi,      diff(E1(c1,xi),c1,xi)=E1c1xi, diff(E1(c1,xi),c1)=E1c1,      diff(E1(c1,xi),xi)=E1xi, E1(c1,xi)=E1,      diff(E2(c2,xi),c2,c2)=E2c2c2, diff(E2(c2,xi),xi,xi)=E2xixi,      diff(E2(c2,xi),c2,xi)=E2c2xi, diff(E2(c2,xi),c2)=E2c2,      diff(E2(c2,xi),xi)=E2xi, E2(c2,xi)=E2];

To save memory we shall use these substitutions and back substitutions for the exponential functions. The substitutions enable us to compute all exponential functions only once in one iteration step.

>	Esu:=[exp(xi/c1)=E1(c1,xi),exp(xi/c2)=E2(c2,xi)]; EBS:=subs(DE,map(u->rhs(u)=lhs(u),Esu));

>	Ebs:=subs(DE,subs(DE,Esu),simplify(remove(has,map(u->rhs(u)=lhs(u),      [diff(Esu,c1,c1)[],diff(Esu,c2,c2)[],       diff(Esu,xi,xi)[],       diff(Esu,c1,xi)[],diff(Esu,c2,xi)[],       diff(Esu,c1)[],diff(Esu,c2)[],       diff(Esu,xi)[]]),0)));

>	s:=subs(Esu,s1+s2);  The whole sum of the residuals - for both functions

>	sm:=selectfun(s,sum);  There are 18 summations, but only 10 of them are unique. We shall collect them

>	s0:=remove(has,sm,[c1,c2]);  We shall compute coefficients independent summations separately

>	Su0:=[seq(s0[j]=Z||j,j=1..nops(s0))]; Bs0:=map(u->rhs(u)=lhs(u),Su0);  Substitutions and back substitutions for the coefficient independent sums

>	s1:=select(has,sm,c1); s2:=select(has,sm,c2);  The coefficient dependent summations

>

DU:=[diff(U1(c1),c1,c1)=U1c1c1,diff(U2(c1),c1,c1)=U2c1c1,      diff(U3(c1),c1,c1)=U3c1c1,diff(U1(c1),c1)=U1c1,      diff(U2(c1),c1)=U2c1,diff(U3(c1),c1)=U3c1,      U1(c1)=U1,U2(c1)=U2,U3(c1)=U3]; DV:=[diff(V1(c2),c2,c2)=V1c2c2,diff(V2(c2),c2,c2)=V2c2c2,      diff(V3(c2),c2,c2)=V3c2c2,diff(V1(c2),c2)=V1c2,      diff(V2(c2),c2)=V2c2,diff(V3(c2),c2)=V3c2,      V1(c2)=V1,V2(c2)=V2,V3(c2)=V3];  Abbreviations for the summations and their derivatives. During iteration each parameter will be assigned to a numerical value, so we have to substitute derivatives by simple names.

>	Su1:=[seq(s1[j]=U||j(c1),j=1..nops(s1))]; Su2:=[seq(s2[j]=V||j(c2),j=1..nops(s2))]; Substitutions for the coeffcient dependent sums

>	Bs1:=subs(DU,map(u->rhs(u)=simplify(lhs(u)),        [diff(Su1,c1,c1)[],diff(Su1,c1)[],Su1[]])); Bs2:=subs(DV,map(u->rhs(u)=simplify(lhs(u)),        [diff(Su2,c2,c2)[],diff(Su2,c2)[],Su2[]]));  and their back substitutions

>	S:=subs(Su0,Su1,Su2,s);  Substituted sum of the squared residuals

>

DA:=[diff(A1(c1,c2,xi),c1)=A1c1,diff(A1(c1,c2,xi),c2)=A1c2,      diff(A1(c1,c2,xi),xi)=A1xi,A1(c1,c2,xi)=A1,      diff(A2(c1,c2,xi),c1)=A2c1,diff(A2(c1,c2,xi),c2)=A2c2,      diff(A2(c1,c2,xi),xi)=A2xi,A2(c1,c2,xi)=A2];  During iteration we have to compute derivatives of the remaining two linear parameters  a1  and a2 . So we shall introduce substitutions for them and for their derivatives. We have to remember that these linear parameters will be a function of the nonlinear parameters c1 , c2  and .

>	a12:=subs(a1=A1(c1,c2,xi),a2=A2(c1,c2,xi),        solve({diff(S,a1),diff(S,a2)},{a1,a2}));  The remaining linear parametrs a1  and a2  are computed from the corresponding normal equations.

>	A12:=subs(DE,DU,DV,DA,a12);  If all preceeding substitutions are used, the result looks simpler.

>	A12c1:=subs(DE,DU,DV,DA,diff(a12,c1)):

>	A12c2:=subs(DE,DU,DV,DA,diff(a12,c2)):

>	A12xi:=subs(DE,DU,DV,DA,diff(a12,xi)):  These derivatives will be used to create a set of three equations needed for Newton's method. Outputs are too long, so they are not displayed.

>

>	Sc1:=subs(a1=A1(c1,c2,xi),a2=A2(c1,c2,xi),diff(S,c1)):  The normal equation  =  for c1 . This equation is extremely nonlinear, so we shall use its linearized shape to compute c1, c2  and  fitting to this equation and similar equations  and  numerically, using Newton's method.

>	SC1C1:=subs(DE,DU,DV,DA,diff(Sc1,c1)):

>	SC1C2:=subs(DE,DU,DV,DA,diff(Sc1,c2)):

>	SC1XI:=subs(DE,DU,DV,DA,diff(Sc1,xi)):

>	SC1:=subs(DE,DU,DV,DA,Sc1):

>

>	Sc2:=subs(a1=A1(c1,c2,xi),a2=A2(c1,c2,xi),diff(S,c2)): The normal equation  =  for c2 .

>	SC2C1:=subs(DE,DU,DV,DA,diff(Sc2,c1)):

>	SC2C2:=subs(DE,DU,DV,DA,diff(Sc2,c2)):

>	SC2XI:=subs(DE,DU,DV,DA,diff(Sc2,xi)):

>	SC2:=subs(DE,DU,DV,DA,Sc2):

>

>	Sxi:=subs(a1=A1(c1,c2,xi),a2=A2(c1,c2,xi),diff(S,xi)): The normal equation  =  for .

>	SXIC1:=subs(DE,DU,DV,DA,diff(Sxi,c1)):

>	SXIC2:=subs(DE,DU,DV,DA,diff(Sxi,c2)):

>	SXIXI:=subs(DE,DU,DV,DA,diff(Sxi,xi)):

>	SXI:=subs(DE,DU,DV,DA,Sxi):

>

>

>

>

>	err:=evalf(rand(-500..500)/25.): Data preparation

>	T:=[seq(j*2.,j=0..200),seq(j*5.,j=81..200)]: N:=nops(T);

>	M:=[seq(evalf(subs(t=T[j],g+err()),4),j=1..N)]:

>	Digits:=20:

>	P1:=plot([seq([T[i],M[i]],i=1..N)],style=point,symbol=circle,color=gray):

>	P1;  Plot of given data

>	c1:=250.: c2:=450.: xi:=470.: dc1i:=1.; dc2i:=1.; dxii:=1.; n:=0; Start:=time();  the initial values for the iteration

>	while max(abs(dc1i/c1),abs(dc2i/c2),abs(dxii/xi))>10^(-5) do;  Iteration procedure   n:=n+1:  Count of the steps   T1:=map(u->`if`(u<=xi,u,NULL),T):  Splitting of data into data sets for the first and second function

>	N1:=nops(T1):

>	M1:=M[1..N1]:

>	T2:=T[N1+1..N]:

>	M2:=M[N1+1..N]:

>	N2:=nops(T2):   EBF:=evalf(EBS);

>	Ebf:=[subs(EBF,Ebs)[],EBF[]];  Substitution of the numerical values of exponential functions and summations and their derivatives     Bf0:=evalf(Bs0);   Bf1:=evalf(Bs1);   Bf2:=evalf(Bs2);

>	A12s:=subs(Ebf,Bf0,Bf1,Bf2,A12):  The numerical values of a1  and a2 and their derivatives

>	A12c1s:=subs(Ebf,Bf0,Bf1,Bf2,A12c1):

>	A12c2s:=subs(Ebf,Bf0,Bf1,Bf2,A12c2):

>	A12xis:=subs(Ebf,Bf0,Bf1,Bf2,A12xi):

>	SC1s:=subs(Ebf,Bf0,Bf1,Bf2,A12s,SC1):

>	SC1C1s:=subs(Ebf,Bf0,Bf1,Bf2,A12s,A12c1s,SC1C1):

>	SC1C2s:=subs(Ebf,Bf0,Bf1,Bf2,A12s,A12c2s,SC1C2):

>	SC1XIs:=subs(Ebf,Bf0,Bf1,Bf2,A12s,A12xis,SC1XI):

>	EQ1:=SC1s+SC1C1s*dc1+SC1C2s*dc2+SC1XIs*dxi:  linearized, dc1 , dc2  and dxi  are corrections, (increments) to the c1, c2 and .

>	SC2s:=subs(Ebf,Bf0,Bf1,Bf2,A12s,SC2):

>	SC2C1s:=subs(Ebf,Bf0,Bf1,Bf2,A12s,A12c1s,SC2C1):

>	SC2C2s:=subs(Ebf,Bf0,Bf1,Bf2,A12s,A12c2s,SC2C2):

>	SC2XIs:=subs(Ebf,Bf0,Bf1,Bf2,A12s,A12xis,SC2XI):

>	EQ2:=SC2s+SC2C1s*dc1+SC2C2s*dc2+SC2XIs*dxi:  linearized, dc1 , dc2  and dxi  are corrections, (increments) to the c1, c2 and .

>	SXIs:=subs(Ebf,Bf0,Bf1,Bf2,A12s,SXI):

>	SXIC1s:=subs(Ebf,Bf0,Bf1,Bf2,A12s,A12c1s,SXIC1):

>	SXIC2s:=subs(Ebf,Bf0,Bf1,Bf2,A12s,A12c2s,SXIC2):

>	SXIXIs:=subs(Ebf,Bf0,Bf1,Bf2,A12s,A12xis,SXIXI):

>	EQ3:=SXIs+SXIC1s*dc1+SXIC2s*dc2+SXIXIs*dxi:  linearized, dc1 , dc2  and dxi  are corrections, (increments) to the c1, c2 and .

>	sol:=solve({EQ1,EQ2,EQ3},{dc1,dc2,dxi}):  The corrections for c1 , c2  and

>	xi:=xi+subs(sol,dxi); dxii:=subs(sol,dxi):

>	c1:=c1+subs(sol,dc1); dc1i:=subs(sol,dc1):

>	c2:=c2+subs(sol,dc2); dc2i:=subs(sol,dc2):

>	FT1:=subs(a1=A1,a2=A2,A12s,f1(t)):  Prepare functions for a plot of the iteration step     FT2:=subs(a1=A1,a2=A2,A12s,f2(t)):

>	Y0:=evalf(subs(t=xi,FT1)):

>	Title:=cat(`Step=`,n,`,  Time=`,convert(time()-Start,string),`sec`):

>	display({plot(FT1,t=0..xi,color=red,thickness=3),            plot(FT2,t=xi..T[N],color=blue,thickness=3),            plot([[0,Y0],[xi,Y0],[xi,0]],color=black),P1},title=Title):   print(%);

>	print(evalf(['c1'=c1,'dc1'=dc1i,'c2'=c2,'dc2i'=dc2i,'xi'=xi,'dxii'=dxii][],8));

>	print(`Max. relative step`= max(abs(dc1i/c1),abs(dc2i/c2),abs(dxii/xi))); od:  If the relative corrections are too big, repeat iteration

>	alias(H=Heaviside):  Abbreviation for the Heaviside  function

>	FT:=FT1*H(xi-t)+FT2*H(t-xi);  The resulting function

>	plot(FT-g,t=T[1]..T[N],title=`Mass difference`);  Difference between new function FT  and Gompertz function g

>	plot({diff(FT,t),diff(g,t)},t=T[1]..T[N],title=`Mass increment`);

>	Mf:=evalf(map(u->subs(t=u,FT),T),4):  Numerical values of the new function FT  at times

>	with(stats):

>	Correlation:=describe[linearcorrelation](Mf,M);  Coefficient of the linear correlation between measured mass  and computed mass .

>

Acknowledgement:

Author would like to thank to Autocont Brno and Waterloo Maple for their continuous interest and support.