3
$\begingroup$

Let's say I have some experimental data (x,y) and I want to find some fitting parameters. But in my mathematical model, the fitting parameter are not in a simple equation that relates x and y. For finding y we should solve some sets of algebraic and differential equations and fitting parameters appears in one of these equations. So currently I am solving these equations using an incremental-iterative numerical method to obtain y. how can I do a curve fit while all the methods in Mathematica want a explicit format of a function.

Edit: The problem I am trying to solve is the growth of damage in a material under cyclic loading. The equations that should be solved are following: enter image description here

where bold variables are 3*3 tensors. These equations show a plasticity model coupled with damage. The numerical method used to solve this initial value problem is called elastic predictor-plastic corrector. The curve fitting parameters that I am trying to find are s and S in the last equation. The discretized from of these equations is as follows: enter image description here

First some material properties are introduced:

Remove["Global`*"]
(* Material Properties *)
EYoung = 10000;
ν = 0.3;
smat = 8.67207899548;
S = 0.0442223584949;
σmeso = 42.4521112168564;
Dc = 0.5;
β = 0.476190476190476;
G = EYoung/(2 (1 + ν));
Cy = 100;
σy = 30;
σf = 10;
h = 0.2;

Two fourth order tensors that show the stiffness and compliance of material are formed:

DElastic = 
Table[  EYoung/(
2 (1 + ν))* (KroneckerDelta[i, l]*KroneckerDelta[j, k] + 
   KroneckerDelta[i, k]*KroneckerDelta[j, l] ) + (
 EYoung*ν)/((1 + ν) (1 - 2 ν)) KroneckerDelta[i, j]*
 KroneckerDelta[k, l], {i, 1, 3}, {j, 1, 3}, {k, 1, 3}, {l, 1, 3}];

CElastic = 
Table[(1 + ν)/(
 2 EYoung) (KroneckerDelta[i, l]*KroneckerDelta[j, k] + 
   KroneckerDelta[i, k]*KroneckerDelta[j, l] ) - ν/
 EYoung KroneckerDelta[i, j]*KroneckerDelta[k, l], {i, 1, 3}, {j, 
1, 3}, {k, 1, 3}, {l, 1, 3}];

Loading on material is cyclic. One cycle of this loading is enough for curve fitting problem:

LoadingHistory = Table[ 1/2*(Sin[x + 3 π/2] + 1), {x, 0, 2 π, 2 π/80.0}];

Three vectors save the evolution of some variables during each increment of loading. The last one is damage which is what we are looking for:

ϵμp =   Table[0, {i, 1, 3}, {j, 1, 3}, {k, 1, Length[LoadingHistory]}];
Xμ = Table[0, {i, 1, 3}, {j, 1, 3}, {k, 1, Length[LoadingHistory]}];
Damage = Table[0, {k, 1, Length[LoadingHistory]}];

Some function are defined to be used:

Devaitoric[tensor_] := Module[{H, tensorD}, H = 1.0/3 Tr[tensor];
tensorD = tensor - 1.0/3 Tr[tensor] IdentityMatrix[3]];

MisesNorm[tensor_] := (
tensorD = Devaitoric[tensor];
Sqrt[3.0/2 Total[tensorD*tensorD, 2]]  );
macaulay[x_] := (Abs[x] + x)/2.0;

The main function, solves the set of equations for one step. t1 is the current step where we have the values and t2 is the next step where we are looking to find the values.

mainfun[σmeso_, smat_, S_, t1_, t2_] := (

σ = LoadingHistory[[t2]]*( {
{0.0, 0, 0},
{0, σmeso, 0},
{0, 0, 0}} );
(* meso strain in integration point -tensor inner product*)
ϵ = 
Table[Sum[
CElastic[[i, j, k, l]]*σ[[k, l]], {k, 1, 3}, {l, 1, 
3}], {i, 1, 3}, {j, 1, 3}];

(*Total Strain at micro level *)
ϵμ = ϵ + β*ϵμp[[All, All, 
t1]]; 
ϵμe = ϵμ - ϵμp[[All, All, 
t1]];
σμ = 
Table[Sum[
DElastic[[i, j, k, l]]*ϵμe[[k, l]], {k, 1, 3}, {l, 
1, 3}], {i, 1, 3}, {j, 1, 3}];
s = σμ - Xμ[[All, All, t1]];
seq = MisesNorm[s];

Δp = 0.0;
f = seq - σf;

If[f > -0.001,

\[ScriptCapitalG] = 3 G (1 - β) + Cy (1 - Damage[[t1]]);
Cs = Table[1, {i, 3}, {j, 3}];
Cp = 1;

While[  Total[Cs, 2] > 1.0*^-13 ∨ Abs[Cp] > 1.0*^-13,


mμ = 3/(2 MisesNorm[s])*Devaitoric[s];
Rs = 
s + 2/3 \[ScriptCapitalG] mμ Δp - 
Table[Sum[
DElastic[[i, j, k, l]]*ϵ[[k, l]], {k, 1, 3}, {l, 1, 
3}], {i, 1, 3}, {j, 1, 3}] + 
2 G (1 - β) ϵμp[[All, All, t1]] + 
Xμ[[All, All, t1]];

Rp = MisesNorm[s] - σf;
Cp = (Rp - Total[ mμ* Rs, 2])/\[ScriptCapitalG];
Cs = 2/3 ( Total[ mμ* Rs, 2] - Rp) mμ - (Rs *MisesNorm[s] + 2/3 \ [ScriptCapitalG]*Δp*Total[mμ* Rs, 2] mμ)/(MisesNorm[s] + \[ScriptCapitalG] Δp);

s = s + Cs;
Δp = Δp + Cp;
];

(*Updating variables to be used in next increment*)

ϵμp[[All, All, t2]] = ϵμp[[All, All, t1]] + Δp mμ;
Xμ[[All, All, t2]] =2/3 Cy (1 - Damage[[t1]]) ϵμp[[All, All, t2]];
σμ = s + Xμ[[All, All, t2]];
σplus = Select[Eigenvalues[σμ], # >= 0 &];
σminus = Select[Eigenvalues[σμ], # < 0 &];
Yμ = (1 + ν)/(2 EYoung) (Total[ σplus*σplus] + h ( (1 - Damage[[t1]])/(1 - h Damage[[t1]]))^2 Total[ \σminus*σminus]) - ν/(2 EYoung) (macaulay[ Tr[ σμ]]^2 +h ( (1 - Damage[[t1]])/(1 - h Damage[[t1]]))^2 macaulay[ -Tr[ σμ]]^2);

Damage[[t2]] = Damage[[t1]] + (Yμ/S )^smat Δp;
, 
Damage[[t2]] = Damage[[t1]];
ϵμp[[All, All, t2]] = ϵμp[[All, All, t1]];
Xμ[[All, All, t2]] = Xμ[[All, All, t1]];
];
Damage
);

Applying this function to the whole loading history will give us the solution over the whole loading:

Do[mainfun[σmeso, smat, S, i, i + 1], {i, 1, 
Length[LoadingHistory] - 1}]

Damage evolution and final value of damage:

ListPlot[Damage]
Last@Damage

The experimental data that we are trying to fit this model to are as follows. DExperiment is the damage after one cycle of loading, which we have tried to find above:

σmesovalues = {32.39618739811888300000, 
32.88596552115932800000, 33.41106740534313200000, 
33.97317756735912300000, 34.57921776339298200000, 
35.23385123371712500000, 35.95417637218495300000, 
36.74066188224332300000, 37.60823476195416500000, 
38.57504805727071300000, 39.66470304378003000000, 
40.96920524497549800000, 42.45211121685638000000};
DExperiment = {1.99976*10^-6, 4.08354*10^-6, 7.67889*10^-6, 
0.00001303, 0.0000212061, 0.0000310106, 0.0000452827, 0.0000633068,
0.0000858743, 0.000110066, 0.000135709, 0.000169626, 0.000196569};
$\endgroup$
  • $\begingroup$ How about using an interpolation function with the points forming the interpolation function being the fitting parameters? $\endgroup$ – Hugh Jan 19 '15 at 6:51
  • 1
    $\begingroup$ What about giving a simple example with the analogous behavior? $\endgroup$ – Alexei Boulbitch Jan 19 '15 at 7:51
  • 1
    $\begingroup$ Test data suggestion from this open question: Not Fitting well exponent BY Mathematica $\endgroup$ – Chris Degnen Jan 19 '15 at 11:18
  • 1
    $\begingroup$ The code that I am dealing with is quite big with many equations. I am trying to trim it down to a simpler code. I will post the results as soon as possible. $\endgroup$ – Ahmad A Jan 19 '15 at 11:50
  • 1
    $\begingroup$ Related question (practically speaking, not conceptually): (47539) $\endgroup$ – Oleksandr R. Jan 19 '15 at 12:46
10
$\begingroup$

I am going to fit an interpolation function to the data. The interpolation points can be adjusted to make a best fit. First we make some data

SeedRandom[12345];
data = Table[{x, Sin[2 \[Pi] x] + 0.1 RandomReal[]}, {x, 0, 0.8,0.02}];

Now we define the $x$ locations which we will use for the interpolation function:

nOfControlPoints = 9;
controlPoints = 
  Subdivide[#1, #2, nOfControlPoints - 1] & @@ MinMax[data[[;; , 1]]];

Here is a plot of the data with the control points showed:

ListPlot[data, Epilog -> {Red, Point[{#, 0}] & /@ controlPoints}, 
 Frame -> True, Axes -> False]

Mathematica graphics

Now we construct the model and fit it with NonlinearModelFit:

model[y : {__Real}] := Interpolation[Transpose[{controlPoints, y}]];
nlm = NonlinearModelFit[data, model[Array[y, nOfControlPoints]][x], 
   Array[y, nOfControlPoints], x];

A comparison with original data:

Plot[nlm["Function"][x], {x, 0, 0.8}, Epilog -> {Point[data]}]

plot

$\endgroup$
  • $\begingroup$ Thanks for the improvements. Much cleaner now and some new code for me to understand. $\endgroup$ – Hugh Jun 6 '15 at 7:25
  • $\begingroup$ Where are Subdivide and MinMax defined? (I can guess what they are, but it would be better to have their actual definitions.) $\endgroup$ – kjo Jun 7 '15 at 14:47
  • 1
    $\begingroup$ @kjo Subdivide and MinMax are both new in version 10.1. I do not have them either since I am still on version 10 (I must update). Alexay Popkov made some useful streamlining to my code and used these. If you look up their definition online they are straightforward to implement. If you have difficulty let me know and I will try to help. Sorry this is a bit unhelpful (for both of us). $\endgroup$ – Hugh Jun 7 '15 at 15:01
  • $\begingroup$ No worries; I ended up defining Subdivide[a_, b_, n_] := Append[Prepend[Table[a + i * ((b - a)/n), {i, n - 1}], a], b] and MinMax[lst_] := {Min[lst], Max[lst]}, and it all worked fine (i.e. I was able to replicate the posted results). $\endgroup$ – kjo Jun 7 '15 at 17:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.