1
$\begingroup$

I have the following system:

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

$\endgroup$
2
$\begingroup$

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)}
$\endgroup$

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.