Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm currently working on solving an SEIR model with data from the CDC. I've created a function to calculate the Square difference for the model and want to create a table of values for the 4 variable combinations, along with their Goodness of fit characteristic from the square difference.

I'm getting this alert when running the function:
"NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended."

It sounds like something I could ignore, but I don't want to mess anything up by doing so.

Also... the big question. Is there any recommendations to make creating this table of values more efficient rather than waiting a half hour for this to work itself out?

ClearAll["Global`*"]

CDCdata = {{0, 161}, {1, 197}, {2, 298}, {3, 339}, {4, 409}, {5, 
614}, {6, 848}, {7, 1359}, {8, 1985}, {9, 3346}, {10, 4590}, {11, 
5938}, {12, 6511}, {13, 6514}, {14, 6952}, {15, 5982}, {16, 
5167}, {17, 4077}, {18, 3126}, {19, 2457}, {20, 2136}, {21, 
1848}, {22, 1680}, {23, 1497}, {24, 1269}, {25, 1007}, {26, 
898}, {27, 644}, {28, 459}, {29, 380}};

ListPlot[CDCdata]

p = Interpolation[CDCdata];
Plot[p[t], {t, 0, 29}]

Population = 20328;
Squarediff[Mu_, Beta_, Alpha_, 
  Nu_] := (sol = 
   First@NDSolve[{s'[t] == 
       Mu*Population - Mu*s[t] - Beta*(i[t]/Population)*s[t], 
      e'[t] == Beta*(i[t]/Population)*s[t] - (Mu + Alpha) e[t], 
      i'[t] == Alpha*(e[t]) - (Nu + Mu)*(i[t]), 
      r'[t] == Nu*i[t] - Mu*r[t],
      i[0] == 161,
      s[0] == 20328,
      r[0] == 0, e[0] == 0},
     {s, e, i, r}, {t, 0, 29}]; 
  Sum[(((CDCdata[[j, 2]] - i[CDCdata[[j, 1]]]) /. sol)^2), {j, 1, 
    Length[CDCdata]}])

values = Flatten[
   Table[{{a, b, c, d}, Squarediff[a, b, c, d]}, {a, 0, 0.0005, 
     0.00001}, {b, 0, 0.2, 0.01}, {c, 0, 0.1, 0.01}, {d, 0, 0.1, 
     0.01}], 1];
share|improve this question
1  
What about Beta? –  Artes Nov 10 '13 at 4:21
2  
Don't create symbols that start with a capital letter as those might shadow/conflict with Mathematica defined symbols. In your case Beta is a Mathematica defined function! –  Matariki Nov 10 '13 at 6:19
    
@Artes are you trying to say the same thing as the comment below? –  Kevin Murphy Nov 10 '13 at 14:26
add comment

1 Answer 1

EDIT

ParametricNDSolveValue may be a useful way to explore the parameter space and then you can quantify goodness of fit. Using the same definitions as your original post (with change of variables):

fun = ParametricNDSolveValue[{s'[t] == 
    mu*Population - mu*s[t] - beta*(i[t]/Population)*s[t], 
   e'[t] == beta*(i[t]/Population)*s[t] - (mu + alpha) e[t], 
   i'[t] == alpha*(e[t]) - (nu + mu)*(i[t]), 
   r'[t] == nu*i[t] - mu*r[t], i[0] == 161, s[0] == 20328, r[0] == 0, 
   e[0] == 0}, {s, e, i, r}, {t, 0, 29}, {mu, beta, alpha, nu}]

This creates set of interpolating functions for: susceptible, infected etc. As you are interested in i[t] this is the third function. You can dynamically explore the parameter space ( I have put the funcion p[t] as the reference):

Manipulate[Module[{foi},
  foi = fun[a, b, c, d][[3]];
  Plot[{foi[t], p[t]}, {t, 0, 29}]],
 {{a, 0.001, \[Mu]}, 0.001, 1, Appearance -> "Labeled"},
 {{b, 0.01, \[Beta]}, 0.01, 1, Appearance -> "Labeled"},
 {{c, 0.01, \[Alpha]}, 0.01, 1, Appearance -> "Labeled"},
 {{d, 0.01, \[Nu]}, 0.01, 1, Appearance -> "Labeled"}
 ]

enter image description here

Not rigorous but a starting point. This may narrow the parameter space and you can calculate or use other methods to find best fit.

share|improve this answer
    
Thanks for your answer! I've changed the variable names. The issue I have now is, I want to get all 4 variable names a bit more accurate. By changing the step size in the values table, I increase the calculation time by a ton. It ran over night and still isn't done. Any recommendations? –  Kevin Murphy Nov 10 '13 at 14:25
    
@KevinMurphy see edit –  ubpdqn Nov 11 '13 at 1:52
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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