Ask Your Question
3

plotting autonomous differential equations

asked 2018-09-25 02:00:20 -0600

AndyH gravatar image

updated 2018-09-25 03:18:11 -0600

I'm trying to plot the slope field and some solutions to

$$\frac{dx}{dt}=x^2-4x$$

I was able to get the slope field to plot with:

x,y=var('x','t')
plot_slope_field(x^2-4*x,(t,0,5),(x,-4,8))

I can't figure out how to the solutions to plot. I can get a general solution:

t = var('t')
x = function('x')(t)
f=desolve(diff(x,t) == x^2-4*x,[x,t])
show(f)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{4} \, \log\left(x\left(t\right) - 4\right) - \frac{1}{4} \, \log\left(x\left(t\right)\right) = C + t$$

If I try it with initial values, say ics=[2,2], I get imaginary solutions:

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{4} \, \log\left(x\left(t\right) - 4\right) - \frac{1}{4} \, \log\left(x\left(t\right)\right) = \frac{1}{4} i \, \pi + t - 2$$

I was able to plot solutions with another windows tool, winplot, but really want to do it in sage. Any clues on how I can do that?

edit retag flag offensive close merge delete

2 answers

Sort by ยป oldest newest most voted
1

answered 2018-09-26 11:30:21 -0600

Emmanuel Charpentier gravatar image

updated 2018-09-26 11:38:14 -0600

Let's see (with a slight change of notations) : we start the same way you did :

sage: y=function("y")(t)
sage: DE=diff(y,t)==y^2-4*y
sage: foo=desolve(DE,y);foo
1/4*log(y(t) - 4) - 1/4*log(y(t)) == _C + t

This is only an _implicit_ solution ; not quite what we seek... Now, Sage seems reructant to transform log(a)-log(b) into log(a/b). With some reason. If we are ready to assume that, _for all t,_ y(t)!=0, we can rewrite this a bit. Let's check our substitution first :

sage: bool(foo.lhs()==1/4*(log((y(t)-4)/y(t))))
True

Okay. We can now get another equation, equivalent to the first under our assumption :

sage: bar=1/4*(log((y(t)-4)/y(t)))==foo.rhs();bar
1/4*log((y(t) - 4)/y(t)) == _C + t

Can Sage solve this ?

sage: solve(bar, y(t))
[log(y(t)) == -4*_C - 4*t + log(y(t) - 4)]

Not yet here.But the exponential has some interesting properties... Can we solve exp(bar.lhs())==exp(bar.rhs()) ?

sage: gee=bar.lhs().exp()==bar.rhs().exp();gee
((y(t) - 4)/y(t))^(1/4) == e^(_C + t)
sage: solve(gee,y(t))
[y(t) == -4/(e^(4*_C + 4*t) - 1)]

Now, we have _an_ explicit solution to your original equation. But there remain a couple of problems to solve :

  • Does this describe _all_ solutions to the original equation ?

  • Under which condition(s) are these solutions _real_ ?

  • Are all the points of these solutions valid ?

Those questions, dear AndyH, are probably the crux of what I guess to be your _homework_, and therefore left to you as the point of the exercise...

_Computation may become easy, analysis is still hard_.

Wisdom of nations...

edit flag offensive delete link more

Comments

Emmanuel,

Thank you for the detailed reply. I didn't know about the lhs, rhs functions. Those are helpful. You are correct this was homework. I'm trying to learn more about sage as I take a differential equations class for fun. It pushes my 58 year old brain to its limits sometimes :)

I cut and pasted your statements, but got slightly different results. I've got 8.3 and 8.4beta6 of sage. In both versions, I get False when I run

bool(foo.lhs()==1/4*(log((y(t)-4)/y(t))))

I also got a depreciation warning. Changing the y(t) to y(t=t) got rid of that. I still got False however. I tried with assume(y(t) != 0), but that didn't make any difference.

solve(bar, y(t))

returned your solution. I didn't need to do the exponentiate it. I don't ...(more)

AndyH gravatar imageAndyH ( 2018-09-29 01:06:50 -0600 )edit
0

answered 2018-10-02 07:50:50 -0600

calc314 gravatar image

Here is another option for plotting. The streamline_plot does a nice job of plotting some sample curves along with arrows to give direction. I like to use the odeint solver since it works for systems so nicely. Here I use it by defining a system $ds/dt = 1, dx/dt = x^2-4x$. This gives the appropriate plot. Note that we have to tell Sage that the number 1 should be in the symbolic ring (SR) to get this to work.

var('x s t')
pp=streamline_plot(x^2-4*x,(t,0,5),(x,-4,5))
initialvalues=[-2,2]
for x0 in initialvalues:
    pp+=line(desolve_odeint([SR(1),x^2-4*x],[0,x0],srange(0,5,0.05),[s,x]),color='red',thickness=3)
pp.show(ymin=-4,ymax=5)
edit flag offensive delete link more

Comments

Thanks! I'm not sure what the symbolic ring is about, but that graph is great. So much to learn :) Looks like I can use that as a template for a project in my class. The professor's emphasis for the project is to use the graphs to explain the mathematical ideas to a lay person. In this case its modeling sustainable fishing. The equation is different from my original post, but your example will help me explore futher. Thank you.

AndyH gravatar imageAndyH ( 2018-10-02 22:54:20 -0600 )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: 2018-09-25 02:00:20 -0600

Seen: 84 times

Last updated: Oct 02