# ODE w/seasonal forcing term

I am working my way through a mathematical modeling book, and I'm having trouble replicating one of the results.

The Differential equation is given as:

$\frac{dC}{dt} = \frac{F}{V}(C_{in}-C)$, where

$F=$ flow in/out

$V=$ volume

$C_{in}=$ concentration of incoming pollutants

$C=$ total concentration at time t

I was able to replicate the results when $F$ and $C_{in}$ are constant. However, I'm experiencing difficulties when adapting this model with a seasonally forcing term.

For the seasonal flow, the following are given:

$F(t)=10^{6}(1+6sin(2 \pi t))$

$C_{in}=10^{6}(10+10cos(2 \pi t))$

I am supposed to plug these into the original model and arrive at the following family of solutions:

I haven't been successful at this (I'm trying to do it in matlab as well).

Below is my try. And if someone knows how to do it matlab, then a million thanks.

NDSolve[{c'[
t] == ((10^6*(1 + 6*Sin[2*Pi*t]))/
V)*((10^6*(10 + 10*Cos[2*pi*t])) - c[t]), c[0] == 10^7},
c[t], t]

• In the argument to Cos, use Pi instead of pi. Also NDSolve cannot include any non-numeric values, i.e., V must be given a value. Commented Sep 19, 2015 at 15:35
• What is the time scale used to derive the model? I.e. what are the assumptions for the length of the time interval for which the differential equation was derived? This would help selecting the range of t in NDSolve. Commented Sep 19, 2015 at 15:45
• As I understand it, the ideal range would be 10 years. The volume is 28*10^6. Thanks for pointing out the my mistakes Bob. Still no luck. Commented Sep 19, 2015 at 15:47
• Is V the quantity that parameterizes that family of curves? If so, look up ParametricNDSolve and use V as the parameter. Commented Sep 19, 2015 at 16:20
• Can you share what book you are using? Commented Sep 19, 2015 at 17:44

It looks like I can get something very similar as the family of solutions in the question using the following commands.

F[t_] := 10^6 (1 + 6.0*Sin[2 π t]);
V = 28*10^6; tEnd = 8;
fsols =
Table[Block[{sol},
Cin[t_] := 10^6*k*(1.0 + 1*Cos[2 π t]);
sol = NDSolve[{c'[t] == F[t]/V (Cin[t] - c[t]), c[0] == 10^7.},
c[t], {t, 0, tEnd}];
c[t] /. sol[[1]]
], {k, 10, 20}];

Plot[fsols/10^6, {t, 0, tEnd}, PlotRange -> All, AspectRatio -> 1/2]


Note the definition of Cin inside Table -- it is just to include the scaling factor k.

It would be better to do this with Manipulate and see the effects of the other quantities used in the model and using different numerical methods in NDSolve.

(I did experiment varying the volume, but I got solutions looking closer to the expected ones by varying the concentration.)

Update -- obtaining a graph closer to the one in the question

In order to get a graph closer to the in the question we can parametrize both the initial condition and the right hand side of the equation. Here is an example:

V = 28*10^6; tEnd = 10;
fsols =
Table[Block[{sol},
F[t_] := 10^6*(1 + 6.0*Sin[2 π t]);
Cin[t_] := 10^6 (10 + 10*Cos[2 π t]);
sol =
NDSolve[{c'[t] == 6*F[t]/V (Cin[t] - c[t]), c[0] == k*10^7.},
c[t], {t, 0, tEnd}, Method -> Automatic];
c[t] /. sol[[1]]
], {k, 0, 0.6, 0.05}];


Using Manipulate

Here is a Manipulate version of the code above with several controls for the influx concentration and the flow through the volume:

V = 28*10^6;
Manipulate[
DynamicModule[{fsols, c, F, Cin},
F[t_] := 10^6 (1 + 6*Sin[2 π t]);
Cin[t_] := 10^6*(10 + 10*Cos[2 π t]);
fsols =
Table[Block[{sol},
F[t_] := 10^6*(1 + 6.0*Sin[2 π t]);
Cin[t_] := 10^6 (10 + 10*Cos[2 π t]);
sol =
NDSolve[{c'[t] == m*F[t]/V (Cin[t] - c[t]), c[0] == k*10^7.},
c[t], {t, 0, tEnd}, Method -> Automatic];
c[t] /. sol[[1]]
], {k, kMin, kMax, 0.05}];
Plot[fsols/10^6, {t, 0, tEnd}, PlotRange -> {All, All},
AspectRatio -> 1/2]
], {{m, 6, "RHS factor"}, 0., 15,
0.5}, {{kMin, 0, "min initial condition factor"}, 0, 2,
0.01}, {{kMax, 0.6, "max initial condition factor"}, 0, 2,
0.01}, {{tEnd, 8, "time interval (years)"}, 1, 20, 0.5}]


It will produce this Manipulate interface:

Update with R-Shiny interface

In the question there was a request for Matlab code. How about R? (I know how to make an interactive interface for this in R.)

library(shiny)
library(deSolve)

Pi <- 3.14159265
V <- 28.0 * 10.0^6;

Fr <- function(t) { 10.0^6 * ( 1.0 + 6.0 * sin(2*Pi*t) ) }
Cin <- function(t) { 10.0^6 * ( 10.0 + 10.0 * cos(2*Pi*t) ) }

PFunc <- function( t, y, m ) { list( m[1] * Fr(t) / V * ( Cin(t) - y[1] ) ) }

server <- function(input, output) {

output$solutionPlot <- renderPlot({ if ( input$kmin < input$kmax) { for( k in seq( input$kmin, input$kmax, 0.05 ) ) { yini <- c( y1 = k*10^7. ) ysol <- ode( y = yini, func = PFunc, times = seq( 0, input$tend, 0.01 ), parms = input$m, method = "ode45" ) if ( k == input$kmin )  {
plot( ysol, type = "l", which = "y1", lwd = 2, ylab = "concentration", main = "", ylim = c( 0, 1 * input\$kmax * Cin(0) ) )
} else {
lines( ysol, type = "l", lwd = 2 )
}
}
}
})
}

ui <- fluidPage(
sidebarLayout(
sidebarPanel(
sliderInput("m", "RHS factor:", min = 0, max = 15, step = 0.5, value = 6.0 ),
sliderInput("kmin", "min initial condition factor:", min = 0, max = 2, step = 0.01, value = 0 ),
sliderInput("kmax", "max initial condition factor:", min = 0, max = 2, step = 0.01, value = 0.6 ),
sliderInput("tend", "time interval (years):", min = 1, max = 10, step = 0.5, value = 8 )
),
mainPanel( plotOutput("solutionPlot") )
)
)

shinyApp(ui = ui, server = server)


One can run that code with RStudio or Rscript. The output would look like this:

• Thanks Anton, this is great. Much appreciated. It is the concentration and flow that I need to manipulate, not the volume. Sorry for not being clear. BTW, in which block of code would I insert the manipulate? Commented Sep 19, 2015 at 16:42
• @JohnWayne360 You are welcome! Please see my update with Manipulate code. Commented Sep 19, 2015 at 16:56
• @ Anton Antonov I didn't really see it yesterday, but is there a reason why the solution is backwards? Do you think the book has it wrong? Commented Sep 20, 2015 at 17:59
• It is possible that the book had it wrong. But it is good to answer the following questions. Which quantities were varied? In what manner? What is the ODE solver that the author(s) used? Commented Sep 20, 2015 at 19:42
• Initial concentrations were varied (C[0]). I believe they used ODE45, maybe Maple. Commented Sep 20, 2015 at 20:46