MATH564 Homework 1: Mathematical modeling

This first homework is intended as a review of some basic mathematical techniques (algebra, calculus, differential equations) and introduction to working within Mathematica notebooks.

Topics

Basic plot construction

Exercise 2, p15

Sub-point (a)

funcsa={π (r^2), ((4π)/3) (r^2)};
Plot[funcsa,{r,0,4},Frame->True,FrameLabel->{"r","(Subscript[f, a])(r)"},GridLines->Automatic]

Plot[Map[PowerExpand[Log[#]]&,funcsa] /. Log[r]->x,{x,-1,4},Frame->True,FrameLabel->{"log(r)","log((Subscript[f, a])(r))"},GridLines->Automatic]

LogLogPlot[funcsa,{r,0,4},Frame->True,FrameLabel->{"r","(Subscript[f, a])(r)"},GridLines->Automatic]

Sub-point (b)

funcsb={3 (x^5), 5 (x^3)};
Plot[funcsb,{x,0,4},Frame->True,FrameLabel->{"x","(Subscript[f, b])(x)"},GridLines->Automatic]

Exercise 3, p16

funcs={3 (x^2)+5x+7, (2^x),((3 (x^2)+5x+7)/(2^x))};
Plot[funcs,{x,0,4},Frame->True,FrameLabel->{"x","f(x)"},GridLines->Automatic]

Celsius to Fahrenheit (univariate linear dependence)

p=Plot[9c/5+32,{c,0,100},Frame->True,FrameLabel->{"TempC","TempF"},GridLines->Automatic]

Transforming power law to linear dependence through logarithms

The following assigns to the variable y the power law A (x^k)

y=A (x^k)

In contrast, this defines a function f(x)=A (x^k)

f[x_] = A (x^k)

The function can subsequently be evaluated

f[1]

A linear dependence can be obtained through logarithms

PowerExpand[Log[y]]

Note the use of a function to carry out the expansion. Other types of expansions are also available

Expand[(((a+b))^3)]

TrigExpand[Cos[α+β]]

Linear regression

Measured data

Import data on AIDS cases

<<"http://mitran-lab.amath.unc.edu/courses/MATH564/lessons/Table222";Map[Length[#]&,{Time,AIDS,CAC}]

TableForm[{Time,AIDS,CAC},TableDirections->Row,TableHeadings->{{"Year","Cases","Cummulative"}}]

ListPlot[Transpose[{Time,AIDS}],Frame->True,FrameLabel->{"Year","Nr. AIDS cases"},GridLines->Automatic]

ListPlot[Transpose[{Time,CAC}],Frame->True,FrameLabel->{"Year","Nr. AIDS cases"},GridLines->Automatic]

lT= Log[Time-1980]; lC = Log[CAC];
data=Transpose[{lT  ,lC }];
lcac[x_]=Fit[data,{1,x},x]

tcac[x_]=Fit[data,{1,Sin[x],Cos[x]},x]

Plot[{lcac[x],tcac[x]},{x,0,3},PlotLegends->"Expressions"]

fit=Transpose[{lT,lcac[lT]}];
ListPlot[{data,fit},Frame->True,FrameLabel->{"log(Year-1980)","log(CAC)"},Axes->False,GridLines->Automatic,PlotLabel->"LogLinear fit of Cummulative AIDS Cases" ]

Synthetic data

Fit synthetic data to an exponential model y(x)=A (e^(r x)). Synthetic data is generated to test some algorithm as opposed to carrying processing of measurements. As an example, define an exponential model

y[x_,A_,r_] = A Exp[r x]

Generate data for some choice of A, r

Aexact = 1; rexact = 0.2;
data = Table[{x,y[x,Aexact,rexact]},{x,-0.9,0.9,0.2}]

Apply logarithms to obtain a model linear in parameters r, log(A)

logy[x_,r_,logA_]=PowerExpand[Log[y[x,A,r]]] /. Log[A]->logA

Apply the transform to the data

logdata =  Table[{x,Log[y[x,Aexact,rexact]]},{x,-0.9,0.9,0.2}]

Carry out a linear fit to the data to find logA, r

yfit[x_]=Fit[logdata,{1,x},x]

fitcoefs = CoefficientList[yfit[x],x]

logAfit = fitcoefs[[1]]; Afit= Exp[logAfit]; rfit=fitcoefs[[2]];
{Afit,rfit}

Note from the above that the exact coefficients are recovered. Now randomly perturb the data to mimic an experimental measurement error.

ampP = 0.01; cenP=0.5;
logdataP =  Table[{x,Log[y[x,Aexact,rexact]]+ ampP(Random[]-cenP) },{x,-0.9,0.9,0.2}]

Repeat the above linear fit steps for the perturbed data

yfitP[x_]=Fit[logdataP,{1,x},x];
fitcoefsP = CoefficientList[yfitP[x],x];
logAfitP = fitcoefsP[[1]]; AfitP= Exp[logAfitP]; rfitP=fitcoefsP[[2]];
{AfitP,rfitP}

Note that in the presence of a centered random perturbation of amplitude 0.01, values close to the exact model are obtained. Investigate what happens with a larger perturbation, or measurements that are not centered (systematic measurement error)

ampP = 0.1; cenP=0.5;
logdataP =  Table[{x,Log[y[x,Aexact,rexact]]+ ampP(Random[]-cenP) },{x,-0.9,0.9,0.2}];
yfitP[x_]=Fit[logdataP,{1,x},x];
fitcoefsP = CoefficientList[yfitP[x],x];
logAfitP = fitcoefsP[[1]]; AfitP= Exp[logAfitP]; rfitP=fitcoefsP[[2]];
{AfitP,rfitP}

ampP = 0.01; cenP=0.;
logdataP =  Table[{x,Log[y[x,Aexact,rexact]]+ ampP(Random[]-cenP) },{x,-0.9,0.9,0.2}];
yfitP[x_]=Fit[logdataP,{1,x},x];
fitcoefsP = CoefficientList[yfitP[x],x];
logAfitP = fitcoefsP[[1]]; AfitP= Exp[logAfitP]; rfitP=fitcoefsP[[2]];
{AfitP,rfitP}

Classroom task: present plots of the fitted coefficients for random perturbation amplitudes (2^k), k=-3,3, and cenP=0,0.5,1.

Multilinear regression

Two independent variables, one dependent variable

Consider a model of jumping ability as a function of muscle size and tone

c={58.5, 50.0, 59.5, 58.0, 60.5, 57.5, 49.3};
lp={220, 150, 165, 270, 200, 250, 210};
vj={19.5,18.0,22.0,19.0,21.0, 22.0, 29.5};
n=Length[c]

Organize the data into a two-column table and carry out a fit to a linear model with one independent variable

data=Transpose[{c,vj}];
Fit[data,{1,x},x]

A different way to organize the data is through the table function, obtaining the same result

data=Table[{c[[i]],vj[[i]]},{i,1,n}];
Fit[data,{1,x},x]

Organize the data into a three column table with two independent variables, and carry out a multilinear fit

data=Table[{c[[i]],lp[[i]],vj[[i]]},{i,1,n}];
f[x_,y_]=Fit[data,{1,x,y},{x,y}]

The computed fit (least squares approximation) can be used for predictions

f[60.5,200]

A visual representation of the quality of the fit. First construct the plane corresponding to the fit, add a representation of the data points, and then show both.

fitplt=Plot3D[f[x,y],{x,49,60},{y,130,250}];
pntplt = ListPointPlot3D[data];
Show[{pntplt,fitplt}]

Lotka-Volterra system

PreyEq=y'[t]==a y[t]-p y[t] z[t];
PredEq=z'[t]==-b y[t]+q y[t] z[t];
LVparams={a->0.7,b->1,p->1.3,q->1};
ICs={y[0]==1,z[0]==0.5};
LVsystem=Flatten[{PreyEq,PredEq,ICs}/. LVparams]

sol=NDSolve[LVsystem,{y[t],z[t]},{t,0,1.5}][[1]]

Plot[Evaluate[{y[t],z[t]}/. sol],{t,0,1.5},PlotLegends->{"Prey","Pred"}]

ParametricPlot[Evaluate[{y[t],z[t]}/. sol],{t,0,1.5}]

Problems

One course credit point is awarded for correct solution of each problem. Exercises are from the textbook, Mathematical Biology: An Introduction with Maple and Matlab by Ronald Shonkwiler and James Herod

Problem 1: Linear regression

Exercise 2, p26

Exercise 4, p27

Problem 2: Multiple regression

Exercise 3, p35

Problem 3: Differential equations

Exercise 1, p45

Problem 4: Difference equations

Chaos example, p52

Carry out a numerical study of the logistic recurrence relation and display the different qualitative behaviors: unique asymptotic value, two values, chaos.