Ask Your Question
0

desolve (how can i enable numeric?)

asked 2019-03-29 11:25:04 +0100

Martin Mårtensson gravatar image

Hi, I am sorry for asking this, but I have been looking all over the internet now and I can't find any solution to this.

When i am using desolve i sometimes get some crazy outputs which is not relevant for me at all.

For instance i am trying to desolve this equation

y' = 0.00054 \cdot y \cdot (500 - y)

So i white

k = var('k')
x = var('x')
y = function('y')(x)

Then i define my de

de = diff(y,x) == 0.00054*1*y*(500-y)

And i desolve like i have done before with many other equations without any problem

desolve(de,y,[0,90])

But i get this result

-100/27*log(y(x) - 500) + 100/27*log(y(x)) == -100/27*I*pi + x - 100/27*log(410) + 100/27*log(90)

The right result should be:

y = 500/(1+4.556*e^(-0.27*x))

Please help me. I have been readning in the doc section about desolve and many other places, but there is so much stuff i dont understand. I am on the edge to install a subsitute CAS program to help me everytime sagemath fails.

edit retag flag offensive close merge delete

Comments

1

The documentation is not as user friendly as a guide; this is useful. I don't have experience with differential equations but adding contrib_ode=True, show_method=True to desolve tells you it's separable and after adding the lines c = de.variables()[0] sol = solve(de,y) ; sol as suggested in the guide (page 219) gets you closer. One of the experts here will know what to do.

dazedANDconfused gravatar imagedazedANDconfused ( 2019-03-29 16:30:33 +0100 )edit

Thank you so much!! I will look in this guide now :D

Martin Mårtensson gravatar imageMartin Mårtensson ( 2019-03-30 09:38:06 +0100 )edit

1 Answer

Sort by » oldest newest most voted
2

answered 2019-03-30 00:55:21 +0100

Emmanuel Charpentier gravatar image

Well, if you insist on numerical solutions to differential equations, you can always try desolve_rk4 or one of the various solutions documented in the relevant chapter of the documentation. But in this case, it's much easier to solve it first in symbolic terms :

sage: var("a,b,c")
(a, b, c)
sage: y=function("y")
sage: de=diff(y(x),x)==a*y(x)*(b-y(x))
sage: S1=desolve(de,y(x),[0,c], ivar=x, contrib_ode=True);S1
-(log(-b + y(x)) - log(y(x)))/(a*b) == (a*b*x - log(-b + c) + log(c))/(a*b)

This is an implicit solution, that can be easily solved : l''s get rid of constants in the left-hand side, by multimpyong both terms by a*b, then simplify the log terms :

sage: S2=(S1*(a*b)).simplify_log();S2
log(-y(x)/(b - y(x))) == a*b*x + log(-c/(b - c))

Lat's take the exponential of both terms of this equality, then expand logs :

sage: S3=S2.lhs().exp().expand_log()==S2.rhs().exp().expand_log();S3
-y(x)/(b - y(x)) == -c*e^(a*b*x)/(b - c)

It remains a first-degree linear equation, whose (unique) solution is trivial :

sage: S=solve(S3,y(x))[0];S
y(x) == b*c*e^(a*b*x)/(c*e^(a*b*x) + b - c)

This solution can (should)) be checked :

sage: bool(diff(S.rhs(),x)==de.rhs().subs({y(x):S.rhs()}))
True

Substitute your numeric constants :

sage: S.subs([a==0.00054,b==500,c==90])
y(x) == 4500*e^(0.270000000000000*x)/(9*e^(0.270000000000000*x) + 41)

One notes that this solution could be expressed approximately as 500/(e^0.27*x)+5.556) (by dividing both numerator and denominator by 9), not quite the expected result...

How did you come to expect this result ? Wasn't there a copy error ?

edit flag offensive delete link more

Comments

Thank you this is perfect!!

Martin Mårtensson gravatar imageMartin Mårtensson ( 2019-03-30 09:38:32 +0100 )edit

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: 2019-03-29 11:25:04 +0100

Seen: 524 times

Last updated: Mar 30 '19