Ask Your Question
0

How to resolve this problem (differential equations) ?

asked 2016-11-03 10:12:15 -0500

updated 2016-11-05 11:05:59 -0500

Hi ! I wrote a code for a project which concerns differential equations but the problem is I don't have the good curves, even if the code is right (I think). Or it's a problem of syntax but I spend a lot of time searching this, without success. I give you the code accompanied with the formulas and the exit of the code. I can't add hyperlink so I write like this or you can't understand my problem. Sorry for inconvenience !

http://image.noelshack.com/fichiers/2...

http://image.noelshack.com/fichiers/2...

http://image.noelshack.com/fichiers/2...

I'm supposed to get them curves :/ The only problem seems to be black curve (on SageMath) which is the green curve on Excel (third picture) and I don't wanted to put the black curve of Excel (so it's normal if it doesn't appear).

http://image.noelshack.com/fichiers/2...

# y[0] = X (Number of likely ie those not infected) = B - nu*X - c*lmd*X
# y[1] = Y (Number of contagious) = c*X*lmd - Y*(v + nu)
# y[2] = Z (Number of non-infectious HIV) = v*Y*(1 - p) - nu*Z
# y[3] = A (Number of AIDS patients) = p*v*Y - A*(d + nu)
# y[4] = lmd (Proportion of contagious, which made themselves pass HIV on the total population) = (Y*beta)/(X + Y + Z)

# params[0] = nu (Natural death rate) == 0.03125 car nu + d ~= d (1 à 1.33/an)
# params[1] = d (Death rates by disease) == 1 (1/an)
# params[2] = B (Nombre de susceptibles immigrants) == 13333.3/an
# params[3] = v (Number of immigrants may) == 0.2/an
# params[4] = c (Number of sexual partners) == 48 (2 à 6/mois soit 48/an)
# params[5] = p (Proportion of catching HIV) == 0.3 (10 à 30%)
# params[6] = beta (Transmission probability) == 0.0014
# t (Time)

@interact
def interactive_function ( nu = slider (0.001, 0.1, 0.005, default = 0.03125),
                           d = slider (0, 1.5, 0.05, default = 1),
                           B = slider (0, 30000, 500, default = 13333.3),
                           v = slider (0, 1, 0.05, default = 0.2),
                           c = slider (0, 365, 1, default = 48),
                           p = slider (0, 1, 0.05, default = 0.3),
                           beta = slider (0, 1, 0.02, default = 0.014) ) :

  def f_1 (t,y,params) :

    return [ params[2] - params[0]*y[0] - params[4]*y[4]*y[0],
             params[4]*y[0]*y[4] - y[1]*(params[3] + params[0]),
             params[3]*y[1]*(1-params[5]) - params[0]*y[2],
             params[5]*params[3]*y[1] - y[3]*(params[1] + params[0]),
             (y[1]*params[6])/(y[0] + y[1] + y[2]) ]


  T = ode_solver()
  T.function = f_1
  T.algorithm = "rk8pd"
  T.ode_solve (y_0 = [60000,40000,0,0,0.0056],
               t_span = [0,35],
               params = [nu,d,B,v,c,p,beta],
               num_points = 1000)

  f = T.solution
  X = [(x[0],x[1][0]) for x in f]
  Y = [(x[0],x[1][1]) for x in f]
  Z = [(x[0],x[1][2]) for x ...
(more)
edit retag flag offensive close merge delete

Comments

1

To display blocks of code, either indent them with 4 spaces, or select the corresponding lines and click the "code" button (the icon with '101 010').

Can you edit your question to do that?

The way it is now, your code is almost illegible, since - some commented-out code is interpreted as headers, - multiplication signs are interpreted as italics markup, - only the indented parts of your code are displayed as code, - while the other parts are messed up with no line breaks.

slelievre gravatar imageslelievre ( 2016-11-04 12:05:34 -0500 )edit

It's ok, I found the solution anyways ^^ Thanks :)

LordHorus gravatar imageLordHorus ( 2016-11-04 13:59:24 -0500 )edit
1

If you've found one, you can enter an answer to your own question.

The questions and answers on Ask Sage help the person who asked the question.

They also help anyone with a similar question in the future, who can find the answer there.

Having the questions nicely typeset and the answers written out is therefore useful.

slelievre gravatar imageslelievre ( 2016-11-04 14:51:27 -0500 )edit

1 answer

Sort by » oldest newest most voted
0

answered 2016-11-05 01:55:46 -0500

I found the problem, it comes from lmd (in the return of def f_1 (t,y,params)) which isn't a differential equation so I don't had to write this there but declare it alone.

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2016-11-03 10:12:15 -0500

Seen: 83 times

Last updated: Nov 05 '16