Finding a numerical value

Question

I have the following system:

When t < 2.5I know that alfa=0 First I had to write a linear aproximation of the system with S>>k. I did that and I obtained the system:

dX/dt = p * X
dS/dt = (p * X) / y 

I solved the this system in mathematica and I obtained { s -> Function[{t}, C[1]-((-1+E^(p*t)) C[2])/y], x -> Function[{t}, E^(p*t)C[2]] } .

My task is to use the solution for X and the data that I have wich is in the format: {{Time, S, X}, {0, 45.84, 0.15}, {0.15 45.29, 0.19}, {0.3, 44.66, 0.25}} to obtain a numerical value for p and I have a hint: Plot the logarithm of X against the time, and estimate the slope Someone said that the linearized solution is: Log[X[t]]==Log[X[0]]+p*t, the use LinearModelFit to fit my (logged) data to a line, p is the slope.

t1 := 0
s1 := 45.84
x1 := 0.15
t2 := 0.15
s2 := 45.29
x2 := 0.15
t3 := 0.3
s3 := 44.66
x3 := 0.25

Solving the equation for x, for t1=0 gives me the constant value C[2]=0.15 so x[t]==0.15*e^(p*t)

solLin1 = Solve[Log[x2] == Log[x[0]] + p*t2, p]
p->1.57
solLin2 = Solve[Log[x3] == Log[x[0]] + p*t3, p]
p->1.70

My question is if I used correct the function Log[X[t]]==Log[X[0]]+p*t in this way? And I don't understand the part useLinearModelFitto fit my (logged) data to a line,pis the slope

Answer 1

In what follows data contains the list of observations without the header.

Evaluating Part[data, All, {1, -1}] will return a list of pairs {t,X} (time and X value pairs). This list, will be transformed accordingly, in order to be used with LinearModelFit (see below).

The solution for X[t] of the differential system is given by X[t]->Exp[p t]C[2] and can be written as Log[X[t]]==Log[C[2]] + p t (taking logs on each hand side and turning -> into ==).

This last (algebraic) expression can be interpreted as a linear equation in t of the form eg y==a+b t. This form can be estimated using LinearModelFit after applying a simple transformation to the available data.

According to the documentation, LinearModelFit requires the data that are passed to it, to be in the form {{x11, x21, ..., x1n, y1}...} (the dependent variable is the last entry in every row of data).

Since we'll be using the logarithmic transformation of the solution for X[t], we are going to need to transform all the X values (the dependent variable in our regression) into (natural) logs. This is accomplished thus: MapAt[Log, Part[data, All, {1, -1}], {All, -1}].

Finally, evaluating lmf=LinearModelFit[transd,t,t]; lmf["BestFit"], where transd holds the transformed data (see above) produces

-1.90346 + 1.70275 t

The (point estimate of the) slope or $\hat{ρ}$ is equal to 1.703 (approx.).

---

edit

code I've used

Clear[naught]

(* returns the symbols in 'x' appended with '0' eg 'naught[x,y]' evaluates to '{x0, y0}' *)
naught[x__] := Map[ToExpression[StringJoin[ToString[#], "0"]] &, {x}]


ClearAll[lin]

SetAttributes[{lin}, Listable]

(* returns a linear approximation of 'f' around 'naught[x]' *)
lin[f_, x__] := Module[{x0 = naught[x], deriv0, f0, diff0},
  f0 = f @@ x0;
  deriv0 = D[f[x], {{x}}] /. Thread[{x} -> x0];
  diff0 = {x} - x0;
  f0 + Total[diff0 deriv0]
 ]


ClearAll[sep]

SetAttributes[{sep}, Listable]

(* returns a list of coefficients for variables 'x' present in 'f'; 
   the first entry contains constant terms eg 'sep[a+bx+cy]' returns '{a,b,c}' *)
sep[f_, x__] := Module[{coefs = Coefficient[f, {x}]},
  Flatten[{Simplify[f - coefs.{x}], coefs}]
 ]

I have defined the rhs of the differential system as

f1[x_, s_] := p s x/(s + k) - a x
f2[x_, s_] := a (si - s) - p/y s x /(s + k)

Evaluating

lin[{f1, f2}, X, S]

returns the linearised system

{
 (-a + (p S0)/(k + S0)) (X - X0) - a X0 + (p S0 X0)/(k + S0) + (S - S0) (-((p S0 X0)/(k + S0)^2) + (p X0)/(k + S0)), 
 a (-S0 + si) + (S - S0) (-a + (p S0 X0)/((k + S0)^2 y) - (p X0)/((k + S0) y)) - (p S0 (X - X0))/((k + S0) y) - (p S0 X0)/((k + S0) y)
 }

and interpreting the fact that S>>k implies S+k->S along with the fact that since the time index in the provided data is below 2.5, the value of a is inferred to be zero (see question)

q = sep[lin[{f1, f2}, X, S], X, S] /. {k + S0 -> S0, k -> 0, a -> 0} // Simplify

returns

{{0, p, 0}, {0, -(p/y), 0}}

which can be shown to be equal to the derived rhs of the system; indeed

Total[{1, X[t], S[t]} #] & /@ q

returns

{p X[t], -((p X[t])/y)}