Message Boards Message Boards

GROUPS:

Fit data by combination of functions, solutions of differential eq.

Posted 1 month ago
618 Views
|
4 Replies
|
2 Total Likes
|

Hello everyone, I start with Mathematica and I don't know how to fit my data by a linear combination of functions that are solutions of several differential equations. More precisely, my experimental data (named after "data") are typically Absorbance versus time (absorb=f(t)). At each time, this absorbance is the sum of the absorbance of different species which concentration depends on time. For example,

absorb[t] == (((a[t] + ag[t])*0.0003654) + (0.0003272*
      h[t]) + (0.000553*f[t]))*100/a[0]

where a(t), ag(t), h(t) and f(t) designate species concentration. As these species are involved in different chemical equilibria, all the concentrations are interdependent and solution of differential equation representing the rate laws:

Derivative[1][a][t] == -k1*a[t]*g[t] + k2*ag[t], 
Derivative[1][g][t] == -k1*a[t]*g[t] + k2*ag[t],  
Derivative[1][ag][t] == k1*a[t]*g[t] - k2*ag[t] - k3*ag[t], 
Derivative[1][gox][t] == k3*ag[t], 
Derivative[1][h][t] == k3*ag[t] - 2*k4*h[t]*h[t] + 2*k5*f[t]*f[t], 
Derivative[1][f][t] == 2*k4*h[t]*h[t] - 2*k5*f[t]*f[t].

So I need to solve the differential equations, and to find the values of k1 to k5 parameters. (In fact, I can fix some of them to specific values, in order to reduce the number of parameters to be found). Mathematica can't find analytical solutions, therefore I opt for numerical solutions. I try something but it doesn't work:

k1 = 157000
k2 = 70
k3 = 0.0275
Clear[k4, k5]
model = ParametricNDsolveValue[{Derivative[1][a][t] == -k1*a[t]*g[t] +
k2*ag[t], Derivative[1][g][t] == -k1*a[t]*g[t] + k2*ag[t],  
Derivative[1][ag][t] == k1*a[t]*g[t] - k2*ag[t] - k3*ag[t], 
Derivative[1][gox][t] == k3*ag[t], 
Derivative[1][h][t] == k3*ag[t] - 2*k4*h[t]*h[t] + 2*k5*f[t]*f[t], 
Derivative[1][f][t] == 2*k4*h[t]*h[t] - 2*k5*f[t]*f[t], 
absorb[t] == (((a[t] + ag[t])*0.0003654) + (0.0003272*
h[t]) + (0.000553*f[t]))*100/a[0], a[0] == 0.000006, 
g[0] == 0.0038, ag[0] == 0, h[0] == 0, gox[0] == 0, f[0] == 0, 
absorb[0] == 100 *0.0003654}, absorb, {t, 0, 2000}, {k4, k5}]
Fit = FindFit[data, model[k4, k5][t], {k4, k5}, t]

Note that data has been defined before. A part of the data is indicated bellow:

data = {{0.2, 0.0366929}, {0.4, 0.0368321}, {0.6, 0.0368786}, {0.8, 
0.0367551}, {1, 0.036677}, {1.2, 0.036687}, {1.4, 0.0365921}, {1.6,
0.0365438}, {1.8, 0.0364348}, {2, 0.0364835}, {2.2, 
0.036489}, {2.4, 0.0363618}, {2.6, 0.0361545}, {2.8, 
0.0362492}, {3, 0.0362564}, {3.2, 0.0362539}, {3.4, 
0.0361269}, {3.6, 0.0362652}, {3.8, 0.0361286}, {4, 
0.0361022}, {4.2, 0.0360444}, {4.4, 0.0360916}, {4.6, 
0.0360261}, {4.8, 0.0360464}, {5, 0.0359304}, {5.2, 
0.0358969}, {5.4, 0.0359568}, {5.6, 0.036002}, {5.8, 
0.0358563}, {6, 0.0360431}, {6.2, 0.0358403}, {6.4, 
0.0358596}, {6.6, 0.0358963}, {6.8, 0.0359946}, {7, 
0.0358273}, {7.2, 0.0359167}, {7.4, 0.0357477}, {7.6, 
0.0357798}, {7.8, 0.0358253}, {8, 0.0358452}, {8.2, 
0.0358645}, {8.4, 0.035836}, {8.6, 0.0358447}, {8.8, 
0.0357238}, {9, 0.0358621}, {9.2, 0.0358738}, {9.4, 
0.0357445}, {9.6, 0.0356904}, {9.8, 0.0357929}, {10, 
0.0355983}, {10.2, 0.0355699}, {10.4, 0.0357328}, {10.6, 
0.0357016}, {10.8, 0.0357402}, {11, 0.0358079}, {11.2, 
0.0355921}, {11.4, 0.0355836}, {11.6, 0.0356723}, {11.8, 
0.0356181}, {12, 0.0356154}, {12.2, 0.0355639}, {12.4, 
0.0357107}, {12.6, 0.0356061}, {12.8, 0.0355821}, {13, 
0.0356598}, {13.2, 0.0356295}, {13.4, 0.0355047}, {13.6, 
0.0354531}, {13.8, 0.0356424}, {14, 0.0355315}, {14.2, 
0.0353596}, {14.4, 0.0355947}, {14.6, 0.0355897}, {14.8, 
0.0354505}, {15, 0.0356059}, {15.2, 0.03553}, {15.4, 
0.0354226}, {15.6, 0.0354775}, {15.8, 0.0355252}, {16, 
0.0355337}, {16.2, 0.0354188}, {16.4, 0.0354711}, {16.6, 
0.0354321}, {16.8, 0.0354442}, {17, 0.0354882}, {17.2, 
0.0353723}, {17.4, 0.0354554}, {17.6, 0.0354401}, {17.8, 
0.0353133}, {18, 0.0354487}, {18.2, 0.035253}, {18.4, 
0.0352872}, {18.6, 0.0354157}, {18.8, 0.0352881}, {19, 
0.0354663}, {19.2, 0.0353374}, {19.4, 0.0351861}, {19.6, 
0.035292}, {19.8, 0.0353496}, {20, 0.0354485}, {20.2, 
0.0353184}, {20.4, 0.0353444}, {20.6, 0.0352798}, {20.8, 
0.0352445}, {21, 0.0352092}, {21.2, 0.0352332}, {21.4, 
0.0351155}, {21.6, 0.0350688}, {21.8, 0.035244}, {22, 
0.0351599}, {22.2, 0.0352065}, {22.4, 0.0351669}, {22.6, 
0.0350925}, {22.8, 0.0352208}, {23, 0.0351599}, {23.2, 
0.0351167}, {23.4, 0.0352664}, {23.6, 0.0352408}, {23.8, 
0.0351471}, {24, 0.0351475}, {24.2, 0.0352359}, {24.4, 
0.0350241}, {24.6, 0.0350356}, {24.8, 0.0351165}, {25, 
0.0350685}, {25.2, 0.0349941}, {25.4, 0.0351912}, {25.6, 
0.0349834}, {25.8, 0.0350234}, {26, 0.0350935}, {26.2, 
0.0350122}, {26.4, 0.0349314}, {26.6, 0.0350089}, {26.8, 
0.0350968}, {27, 0.0349774}, {27.2, 0.0349673}, {27.4, 
0.0351262}, {27.6, 0.0349908}, {27.8, 0.0350037}, {28, 
0.0350851}, {28.2, 0.0349684}, {28.4, 0.0350994}, {28.6, 
0.0349293}, {28.8, 0.0350352}, {29, 0.0350845}, {29.2, 
0.0349797}, {29.4, 0.0348336}, {29.6, 0.0349395}, {29.8, 
0.0348764}, {30, 0.0349416}, {30.2, 0.034924}, {30.4, 
0.0349941}, {30.6, 0.0349116}, {30.8, 0.0350123}, {31, 
0.0350541}, {31.2, 0.0348924}, {31.4, 0.0348759}, {31.6, 
0.0350295}, {31.8, 0.0348951}, {32, 0.0350272}, {32.2, 
0.0349144}, {32.4, 0.0348181}, {32.6, 0.0350342}, {32.8, 
0.0348983}, {33, 0.0349111}, {33.2, 0.034879}, {33.4, 
0.0349684}, {33.6, 0.0349732}, {33.8, 0.034925}, {34, 
0.0348699}, {34.2, 0.0349309}, {34.4, 0.0349187}, {34.6, 
0.03491}, {34.8, 0.0348181}, {35, 0.0350401}, {35.2, 
0.0349524}, {35.4, 0.03494}, {35.6, 0.0349956}, {35.8, 
0.0349555}, {36, 0.0349272}, {36.2, 0.0349614}, {36.4, 
0.0350246}, {36.6, 0.0349111}, {36.8, 0.0349796}, {37, 
0.0350267}, {37.2, 0.0348849}, {37.4, 0.0349613}, {37.6, 
0.0348261}, {37.8, 0.0347946}, {38, 0.0348967}, {38.2, 
0.0350256}, {38.4, 0.0349309}, {38.6, 0.0348475}, {38.8, 
0.0349582}, {39, 0.0348518}, {39.2, 0.0350048}, {39.4, 
0.0349464}, {39.6, 0.0347781}, {39.8, 0.0349357}, {40, 
0.0349159}, {40.2, 0.0348935}, {40.4, 0.0348758}, {40.6, 
0.0348598}, {40.8, 0.0347577}, {41, 0.0349309}, {41.2, 
0.0348737}, {41.4, 0.034925}, {41.6, 0.034956}, {41.8, 
0.0350867}, {42, 0.0347737}, {42.2, 0.0349256}, {42.4, 
0.0347705}, {42.6, 0.0349143}, {42.8, 0.0347733}, {43, 
0.0349668}, {43.2, 0.0348357}, {43.4, 0.0347946}, {43.6, 
0.0349791}, {43.8, 0.0346116}, {44, 0.0349122}, {44.2, 
0.0349095}, {44.4, 0.0347832}}

I thank you in advance for your advice and all the help you can give me. Sincerely Helene.

4 Replies

I suggest you give a simpler example of this type of problem.

Posted 1 month ago

Here is a simplified version. The heart of the matter remains the same. Experimental data are variations of absorbance versus time which can be described by the following analytical expression:

absorb[t] == (ea*a[t] + eh*h[t] + ef*f[t])*100/a[0]

with ea=0.0003654, eh= 0.0003272 and ef=0.000553 This expression involved 3 species named a, h and f which concentration variation with time is governed by the following differential equations (Mathematica doesn't find any analytical solutions)

NDSolve[{a'[t] == -k3*a[t], h'[t] == k3*a[t] - 2*k4*h[t]*h[t], 
  f'[t] == 2*k4*h[t]*h[t] , a[0] == a0, h[0] == 0, f[0] == 0}, {a, h, 
  f}, t]

To fit my data i tried

k3 = 0.0275
Clear[k4]
sol = ParametricNDSolve[{ Derivative[1][a][t] == -k3*a[t],  
   Derivative[1][h][t] == k3*a[t] - 2*k4*h[t]*h[t], 
   Derivative[1][f][t] == 2*k4*h[t]*h[t], a[0] == 0.000006, 
   h[0] == 0, f[0] == 0}, {a, h, f}, {t, 0, 2000}, {k4}]
model2[k4][
   t] = ((a[t]*0.0003654) + (0.0003272*h[t]) + (0.000553*f[t]))*100/
    a[0];
fit = FindFit[data, {model2[k4][t] /. sol, 0 < k4 < 1000}, {k4}, t];
Plot[Evaluate[model2 /. fit], {t, 0, 850}]

It doesn't work ... I don't well understand how Mathematica can proceed to fit my data. The above code implies that Mathematica has to solve numerically differential equations depending on a parameter (k4) of which it has no idea of the order of magnitude. After it creates the analytical model then it fits the data to find k4. To my opinion, it doesn't make sense. The ideal is that it adjusts and solves the differential equations at the same time to find k4, but I don't know how to write that. Many thanks for your advice Helene

Hi Helene,

I guess basically your definition of model2 is not correct. Try:

k3 = 0.0275;
Clear[k4]
sol = ParametricNDSolve[{
    Derivative[1][a][t] == -k3*a[t],
    Derivative[1][h][t] == k3*a[t] - 2*k4*h[t]*h[t], 
    Derivative[1][f][t] == 2*k4*h[t]*h[t],
    a[0] == 0.000006, h[0] == 0, f[0] == 0}, {a, h, f}, {t, 0, 2000}, {k4}];
{aFunc, hFunc, fFunc} = {a, h, f} /. sol;
model2[k4_][t_] := ((aFunc[k4][t]*0.0003654) + (0.0003272*hFunc[k4][t]) + (0.000553*fFunc[k4][t]))*100/aFunc[k4][0];
fit = FindFit[data, {model2[k4][t], 0 < k4 < 1000}, {k4}, t];
Show[ListPlot[data], Plot[model2[k4 /. fit][t], {t, 0, 50}]]

enter image description here

Hi Henrik, thanks for your answer. It works well as I expect.

Helene

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract