Material parameter estimation using GPL Maxima

Maxima is a [Computer Alegbra System]? which is freely available under the Gnu Public License. Its home page is http://maxima.sourceforge.net .

For many computer algebra tasks, including material law parameter estimation, it is a viable alternative to Mathematica.

There are several frontends available but I recommend wxMaxima http://wxmaxima.sourceforge.net .

Here is my Maxima script for estimating some parameters for the Costa law. It is not optimal but gives an idea of what can be done.

/*
Optimise material parameters for non-linear finite strain elasticity
using Levenberg-Marquardt algorithm.
Glenn Ramsey November 2005. Based on Holger Schmid's Mathematica
implementation.
NOTE Requires maxima > 5.9.2, currently (24/11/05) that means using the cvs version.
*/

kill(all);

/* set up some matrices to represent the deformation gradient tensor */
Fnf: matrix([1,k,0], [0,1,0], [0,0,1])$
Ffn: transpose(Fnf)$
Fns: matrix([1,0,0], [0,1,0], [0,k,1])$
Fsn: transpose(Fns)$
Ffs: matrix([1,0,0], [0,1,0], [k,0,1])$
Fsf: transpose(Ffs)$
Fss: matrix([1,0,0],[0,1,0],[0,0,1+k])$

Nsf: matrix([0,0,1])$
Nsn: matrix([0,0,1])$
Nss: matrix([0,0,1])$

/* F, NA are global lists of matrices */
/* Put the matrices in a list so we can refer to them by index */
F:[Fsf, Fsn, Fss];
NA:[Nsf, Nsn, Nss];

dim:3$
(f:1, n:2, s:3)$
maxstrain:0.25$


/* Use named variables for indices so we know what mode we are dealing with. */
/*
Data is available for 2 shear modes and 1 normal mode so here
I'm setting things up so that the modes we don't have data for are
not used. Doing it this way makes it easier if I change the material
law, then I only have to rearrange the index numbers in one place.
*/
nparam:4;
(ff:nparam+1, fn:nparam+1, fs:1, nn:nparam+1, ns:2, ss:3)$
(aa:4)$

/* Strain energy function - Costa et al, Ph.Tr.R.Soc.Lond.,359,2001 */
/*
Q(b):=b[ff,1]*(E[f,f])^2 + \
2*b[fn,1]*((1/2)*(E[n,f]+E[f,n]))^2 + \
2*b[fs,1]*((1/2)*(E[f,s]+E[s,f]))^2 + \
b[nn,1]*(E[n,n])^2 + \
2*b[ns,1]*((1/2)*(E[n,s]+E[s,n]))^2 + \
b[ss,1]*(E[s,s])^2 ;

W(b):=(1/2)*b[a,1]*(exp(Q(b))-1);
*/

/* The above form doesn't work, so remove unused terms and subst for Q(b) */

W(b):=(1/2)*b[aa,1]*(exp( \
2*b[fs,1]*((1/2)*(E[f,s]+E[s,f]))^2 + \
2*b[ns,1]*((1/2)*(E[n,s]+E[s,n]))^2 + \
b[ss,1]*(E[s,s])^2  )-1);

/* 2nd Piola-Kirchhoff stress */
/* derivative of strain energy function wrt components of strain tensor */
S:zeromatrix(dim,dim)$
for i thru dim do for j thru dim do S[i,j]:diff(W(b),E[i,j])$

/* function for t (model force) in terms of F, N, x */
tm(v,F,N,x):=block([b,k,E],b:v,k:x,Fev:ev(F),E:(1/2)*(transpose(Fev).Fev-ident(dim)),ev(Fev.S.N));

/* Experimental force is given by polynomials fitted to the data
These cubics do not fit the data well for small strains, but this
is what was published...  Douglas et el,JEB,201,1998.
*/

fpd(x):=4970*x^3 + 6320*x^2 - 425*x + 26$
flm(x):=-1090*x^3 + 3800*x^2 - 300*x +21$
fnn(x):=1.39e5*x^3 + 230*x^2 - 522*x + 28$

Te:[fpd,flm,fnn];

te(x,mode):=(Te[mode])(x);

ta(x,mode,v):=tm(v,F[mode],NA[mode],x)[mode,1];

plot2d([(te(x,1)),(te(x,2)),(te(x,3))],[x,0,maxstrain],[y,0,200], \
[gnuplot_term,ps], [gnuplot_out_file, "experimental.ps"]);

/* Generate a list of strain steps for the objective function. */
u:block([y],y:[], for n:0 thru maxstrain step maxstrain/10 do y:append(y,[n]),y)$
uln:length(u)$

/* objective function in terms of material parameters*/

/* display the objective function in noun form */
'sum('sum((ta(u[n],mode,v)-te(u[n],mode))^2,n,1,uln),mode,1,3);

/* now actually define it */
forcediff(mode,v):=sum((ta(u[n],mode,v)-te(u[n],mode))^2,n,1,uln);
obj(v):=sum(forcediff(mode,v),mode,1,3);

/* have to do this substitute and evaluate thing there, not sure why */
m:obj(b)$
mobj(v):=block(b:v,(ev(m)))$

/* partial derivatives of the objective function wrt material parameters */
g:zeromatrix(nparam,1)$
for i:1 thru nparam do g[i,1]:diff(obj(b),b[i,1]);
grad(v):=block(b:v,r:zeromatrix(nparam,1),for i:1 thru nparam do (r[i,1]:ev(g[i,1])),r);

/* 2nd partial derivatives of the objective function wrt material parameters */
h:zeromatrix(nparam,nparam)$
for i:1 thru nparam do for j:1 thru nparam do h[i,j]:diff(g[j,1],b[i,1]);
hess(v):=block(b:v,r:zeromatrix(nparam,nparam),
  for i:1 thru nparam do for j:1 thru nparam do (r[i,j]:ev(h[i,j])),r);

/* Now the equations are set up */

/* Optimization loop Levenberg-Marquardt*/

/* Initial guess */
/* theta[0]:transpose(matrix([20,20,20,20])); */
theta[0]:matrix([17],[9],[17],[90]);
lmda:1;
aom:1;
ath:1;

for i:0 while (aom > 10^-5 or ath > 10^-5) do (
  omega[0]:float(mobj(theta[0])),
  delta:-invert(float(hess(theta[0])+lmda*ident(nparam))).float(grad(theta[0])),
  display(delta),
  theta[1]:theta[0] + delta,
  omega[1]:float(mobj(theta[1])),
  display(omega[1]),
  display(omega[0]),
  display(theta[1]),
  if (omega[1]<omega[0]) then (
    ath:(theta[1]-theta[0]).(theta[1]-theta[0]),
    aom:abs(omega[1]-omega[0]),
    theta[0]:theta[1],
    display(ath),
    display(omega[1]),
    display(theta[1]),
    lmda:lmda/10 ) else (
    lmda:lmda*10 ),
  i:i+1,
  display(i)
);


/*
r:zeromatrix(nparam+1,1);
for i thru nparam do r[i,1]:theta[1][i,1];
*/

r:theta[1];

/* display some results */
te(0.1,1);
float(ta(0.1,1,r));
te(0.1,2);
float(ta(0.1,2,r));
te(0.1,3);
float(ta(0.1,3,r));
plot2d([(te(x,1)),ta(x,1,r),(te(x,2)),ta(x,2,r),(te(x,3)),ta(x,3,r)],[x,0,maxstrain],[y,0,2000], \
[gnuplot_term,ps], [gnuplot_out_file, "fitted.ps"]);