# how to find curve fitting parameters without an explicit relation between x and y?

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: 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: 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}];


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]}];


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

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_] := (

{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
);


Do[mainfun[σmeso, smat, S, i, i + 1], {i, 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};

• How about using an interpolation function with the points forming the interpolation function being the fitting parameters? – Hugh Jan 19 '15 at 6:51
• What about giving a simple example with the analogous behavior? – Alexei Boulbitch Jan 19 '15 at 7:51
• Test data suggestion from this open question: Not Fitting well exponent BY Mathematica – Chris Degnen Jan 19 '15 at 11:18
• 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. – Ahmad A Jan 19 '15 at 11:50
• Related question (practically speaking, not conceptually): (47539) – Oleksandr R. Jan 19 '15 at 12:46

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;
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] 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]}] • Thanks for the improvements. Much cleaner now and some new code for me to understand. – Hugh Jun 6 '15 at 7:25
• Where are Subdivide and MinMax defined? (I can guess what they are, but it would be better to have their actual definitions.) – kjo Jun 7 '15 at 14:47
• @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). – Hugh Jun 7 '15 at 15:01
• 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). – kjo Jun 7 '15 at 17:04