# ODE-system solving with initial conditions (Maple)

This is my problem:

I have to solve a differential-equation-system, and then to plot the results. The system is :

{ S'=-r*S*J,

J'=r*S*J-a*J,

R'=a*J }.

It's about epidemiology, S, J, and R are the Susceptible, Infected, and Retired people, while a and r are two parameters that I determined before.

The thing is, when I give as initial conditions the 3 values at t=0, everything is OK, but I don't know these 3 values ! So I'd like to give Maple other values, as the derivate of one of the functions at some point, or this sort of thing... But Maple doesn't accept, whatever my syntax is. I've tried many ways, and already checked the online help that they provide, but didn't find anything helpful. Does anybody know how to fix that issue?

Thank you :)

• Is R decoupled from the other equations and you can solve for it? – Amzoti Jan 15 '14 at 18:18
• What are the values of the parameters $a$ and $r$ that you are using? And what are the three boundary conditions that you want to use? – Carl Love Jan 15 '14 at 20:47
• @Amzoti : R(t) symbolizes the people who are dead from the disease, and the only thing we know about it is that its speed is proportional to the people who are infected, I(t), by coeficient a. However, thanks to some data I found, I can know R(t) and R'(t) for almost any t. Notice that I don't have a general expression of R, only values from a table. – BusyAnt Jan 15 '14 at 21:19
• @CarlLove : With my data and some calculation I did, I have a=1/3 and r=4.90*10^-6. At best I would use R(0)=0, R(infinity)=40000, and another intermediate value of R, that I know at many points. I also have J(infinity)=0, and S(infinity)=N-R(infinity), with N=90000 the number of people when the epidemy begins. I might find a date where S or J is estimated in the archives. – BusyAnt Jan 15 '14 at 21:27
• Maple will not accept infinity as a boundary for a numeric BVP. We usually get around this by replacing infinity with a large value of the independent variable. So what is a large value of $t$? – Carl Love Jan 15 '14 at 22:02

## 2 Answers

Using $t=100$ as an effective value of $\infty$ and using the other values that you gave in the comment, here's how to code this in Maple:

sys:= {
diff(S(t),t) = -r*S(t)*J(t),
diff(J(t),t) = r*S(t)*J(t) - a*J(t),
diff(R(t),t) = a*J(t),
# Boundary conditions:
S(100) = 50000, R(0) = 0, R(100) = 40000
}:

Sol:= dsolve(eval(sys, [r=4.9e-6, a=1/3]), numeric):

plots:-odeplot(
Sol, [[t,S(t)], [t,J(t)], [t,R(t)]], t= 0..100,
legend= ["susceptible","infected","retired"]
);

• Thank you very much ! I didn't know this syntax, the result is exactly what I was searching for, so you saved me from hours of further research ! – BusyAnt Jan 16 '14 at 12:38
• Can't vote up :/ I don't have enough reputation yet But I accepted the answer now it's done :) – BusyAnt Jan 16 '14 at 21:11
• Boundary-value problems are much trickier to solve numerically than initial-value problems. For example, I can't extend the above to $\infty=200$. But I can use reasonable initial values: $S(0)=90000, J(0)=1, R(0)=0$. This removes the problem of the fake $\infty$. I got reasonable plots with this. – Carl Love Jan 17 '14 at 2:48
• I agree with you on that point : it's far better (I mean easier) to take initial values than boundary ones. But I must model the plague thanks to the data I have. Even if the initial values you gave me make correct plots, I'm not allowed to use them ; it'll be like cheating on the data... Else, I've noticed that the fake infinity couldn't go too high. Even with 150, "Newton iteration is not converging". Do you know any way of fixing this? Thanks again for the interest you give – BusyAnt Jan 18 '14 at 20:41
• The "Newton iteration is not converging" is a problem specific to BVPs, and it is very difficult to fix. There is some limited help available at ?dsolve,numeric_bvp,advanced. Note that this error is different from "initial Newton iteration is not converging." – Carl Love Jan 19 '14 at 6:45

Here's what I tried, following your advice :

sys := {J(0) = 1, R(0) = 0, S(0) = 89999, diff(J(t), t) = r*S(t)*J(t)-a*J(t), diff(R(t), t) = a*J(t), diff(S(t), t) = -r*S(t)*J(t)};

Sol:= dsolve(eval(sys, [r=4.9e-6, a=1/3]), numeric);

plots:-odeplot( Sol, [[t,S(t)], [t,J(t)], [t,R(t)]], t= 50..150, legend= ["Susceptible","Infected","Retired"], color=["Blue", "Red", "Green"]);

As you see, the plot is the same, just shifted to the right for approximately 50 days. And no Newton iteration convergence problem anymore ! I think I'll use this solution : by superimposing it over the previous one, I can easily detremine the "true" initial conditions that I need.

EDIT: I searched for the common maximum for J on the two plots, hence I deduced : OldSol(t)=NewSol(t+40)