Differential equation in Derive and GNUPlot

[Davide82]

Hi!

I am looking for some help using Derive.

I have a differential equation, you can see it in the third image (I wanted to use LaTeX to render it but apparently I can't post LaTeX in this section of the forum).
The value of lambda is given in the first image.

I know the function should be solved by numerical methods since it is complicated.
The function should have a plot like the attached one. But I don't know how I can obtain such a plot with derive. I tried to solve with DSOLVE1 but can't get the plot.

Ideally I would like to render the plot with GNUPLOT.
I would be glad if someone can help me getting that plot with both gnuplot and derive.

Thank you all!

 

[Bill Simpson]

Knowing nothing more about you or what you have managed to get working other than what you have literally written above, it is difficult for me to how much progress you have made, what you have working and where you are stuck.

I don't know whether you are having problems just trying to understand DSOLVE1 or whether you are far beyond that. I don't know which of your variable are actually known numeric constants.

If I start with a much simpler problem and try to get that to work first I choose something like sin(x)-y'=0 with y=1 at x=0. From the DSOLVE1 documentation I simplify DSOLVE1(SIN(x),-1,x,y,0,1) and that gives me COS(x)+y=2. I then SOLVE(COS(x) + y = 2, y) and it gives me y=2-COS(x). I flip over to the 2Dplot window, click Plot and it gives me the resulting plot of 2-cos. I then click on File and Export and see that I can export that as a DIB file, JPEG file, TARGA file or TIFF file. But that is having Derive do the rendering. If you want to have gnuplot do the rendering then you are going to have to find a way scrape the equation off the Algebra window and get it into gnuplot. The details of doing that I can't help you with, I don't claim you can do things that I haven't literally just done at the keyboard.

Now back to your original problem and assume all that was far beneath your understanding. If I try DSOLVE1 with some guesses for tau and s I get ? which probably means I have chosen values where there is no solution.

If I back off to using
DSOLVE1_GEN(x*(255/(t*x^5)*(12 + 6*x + x^2))/(1.13*s)*(e^(-x) - y*(1 + e^(-x))), -1, x, y, c) and Simplify that I get an implicit solution involving an integral of a stack of exponentials. Perhaps with the coefficients for s and tau that can be solved. I can't tell.

If that can't work and you will be satisfied with a numeric approximation then I pull in the ODEApproximation.mth library and consider using TAYLOR_ODE or EULER_ODE or any of the other tools in the library.

But I am getting further and further from knowing what you know, what you have done, what works or where you need to go.

If you can provide substantially more detail about what you have done and exactly what step you are stuck at then perhaps I can take a minute to try to offer a suggestion

 

[Davide82]

Dear Bill,

first of all I thank you for your nice and elaborated answer.
You have been very kind.

I apologize because I realized that my question was very poorly written.

I re-write it here in a clearer way:
My aim is to embed a nice plot in a LaTeX document.
The equation I would like to plot is now more clearly written in the attacched picture...
It comes from the Boltzmann transport equation applied to some particles in the early universe.
Since I read that this equation is a bit complicated and should be approached numerically, I opted for calculating it with the computer, in that way I also learn something along the way...
Of course, if someone knows an analytical solution I would be glad to learn that, too.
The plot I am expecting is exactly the one in the first post.

So, as you suggested, I tried to use DSOLVE1(..., -1, x, y, 1, 0.5) and solve for y...
Then I pushed "Plot" and Derive hangs.
Oddly enough, if I change 0.255 --> 0.25 it draws something without hanging (!) but it is a bit different from the plot I expect. I am puzzled.
For your convenience I attached a Derive DFW file (zipped).

If possible, I would also like to plot it using GNUPlot. I am no expert but I believe GNUPlot can manage simple expressions like y = sin(x) but I think it cannot accept differential equations, so I have to make a point-by-point file and let GNUPlot read it (what is a good way to produce that file?).


Thank you once more!


P.S.: why in this section of the forum we are not allowed to use inline TeX in posts? It is so useful!

 

[Davide82]

Now I made some tests and realized that Derive hangs if I put in the expression any number with 3 decimal digits in place of 0.255. I had to write 0.26...
Then it plots something. I don't know the reason of this behaviour.
But, again, the plot is different from what I expected.
There should be an asymptote for y = 0.151

You can find the equation here: http://nicadd.niu.edu/~bterzic/PHYS652/Lecture_15.pdf

The expected plot can be found here: http://www.ge.infn.it/~gracco/Astr_part_0708/6-Nucl_sintesi-Ricomb.pdf on page 13.

 

[Davide82]

I tried TAYLOR_ODE1 and it plots something.. but at x = 1.25 it increases! No asymptote. I used a 10 order Taylor approximation.
Plot is attached.
Derive file, too.

 

[Bill Simpson]

Ok, we have many issues to deal with here.

1. Since the paper you are reading says this should be attacked numerically, spending time on DSOLVE1 is probably wasted.

2. In your eq.dwf on line 2 you have a d@ and that probably says you have made some error in the syntax of using DSOLVE1.

3. Usually computer algebra packages trying to find exact algebraic solutions will have much greater chance of success if you enter exact rationals and do not use decimal points. The reasons for that are more than I have time to explain.

4. Derive, in my experience, VERY rarely "hangs." Usually long delays are because you have just asked it to calculate something that will take a few minutes or a few hundred hours and it is busy doing exactly what you told it to do. If you can demonstrate an actual bug in Derive I would be very interested, even though this was killed by TI not long after they purchased it. In ten or fifteen years and tens of thousands of hours of intense use I did actually find a few bugs which were promptly fixed before TI got their hands on it.

5. If you are trying to "glue" Derive to Gnuplot I would recommend starting with a very simple example, as I did, use Derive to find either the closed form or an approximate form and then somehow transport that (nondifferential) equation over into Gnuplot. Once you have overcome all the difficulties of gluing these together with a trivial example then you can start trying to do this with your far more complicated actual problem.

6. "Why can't we have Tex, Why can't we desktop publish everything? etc." I really do not understand the overwhelming compulsion that drives people to think they have to turn everything from simple text into a picture to display here, just so that picture has to be manually typed back into the software to try to figure out how to fix your problem. Sorry.

7. .255 versus .26 versus... etc. See item 3 above.

8. "why doesn't my Taylor approximation look like the graph in the journal article?" I presume you realize that Taylor is going to generate a polynomial approximation that will be very good near the point at which you approximate the original equation. I presume you also realize that this polynomial approximation will be worse and worse the farther you get from the original point. Consider approximating Sin(x) by a 10'th order Taylor at 0 and ask why it is a bad approximation at x=20. This is made worse because if you look at the graph in the journal article you can probably see that no reasonable polynomial is going to have that shape, no matter what polynomial you choose.

I hope those clear up a small part of the first layer of things that need to be resolved. Once those have all been resolved then we can begin to attack all the other problems at the first layer and hopefully start in the direction of finding all the problems at the second layer.

 

[Davide82]

Thank you once more for helping me.

I try to address every point in your post so that the discussion has some order:
1. The paper says it is complicated and it is better to solve it numerically; nevertheless if a software could find an analytical solution it would be nice. This is why I tried DSOLVE1.

2. The syntax of DSOLVE1 seems pretty easy; I think I didn't make any mistake using it. The syntax I used is in the first line; I believe it represents correctly the equation I am trying to solve. I don't know why Derive leaves d@.

3. Ok, I don't need that explanation (even if I am interested and you seem to have good knowledge!).

4. I said "hangs" because when asked to plot, it hanged for about one minute and then said "The highlighted expression can't be plot" (or something similar); yes, it is not a rigourous "hang". By the way, my expression doesn't seem that much complicated at least to my untrained eyes, so I hope my PC won't took hours or even minutes to solve it! I am lucky to speak with someone which found some bugs in Derive... you must be a very good user and have a lot of experience!

5. Ok. I will address this issue once I got Derive working with me.

6. I understand your point. Indeed it is better to write equations in plain text to make cut & paste easier. But I was simply worrying about readability of the equation. You, for sure, have no problems, but it is faster to identify an expression if rendered graphically. Since in other sections of the forum users are allowed to use TeX inline I was just wondering why here it is forbidden.

8. Yes I know Taylor tries to approximate around a specific point but I hoped that with a 25 order Taylor approximation I could have a nice approximation a bit further away.


If the Taylor approximation is awful (at decent orders... I tried up to 25 where it tooks some minutes to my PC to complete) and an analytical solution is out of question, how would one get a plot like the one of the first post?
Should we use other numerical methods such as Runge-Kutta? What could I do with Derive?

As you may have understood, I am by no means an expert with differential equations neither with Derive, so forgive me if my questions are very basic.

P.S.: from your words, you seem quite a deep expert about Derive. Few days ago a friend of mine told me that Derive is "mediocre"; he is a Mathematica user. In your opinion what are the strong features of Derive in respect to other pieces of software such as Mathematica, Maple, Matlab or other free alternatives? I know they are very different but why you like Derive a lot?

P.S.2: just out of curiosity... what is your scientific background? or what you do in your life? feel free not to answer! I know I am curious! :D

Thank you very much... I know you are taking me in hand...

 

[Bill Simpson]

2. That @ worries me and that is usually a hint something isn't right. If I look carefully at the @ on your line 2 in eq.dfw I have figured out this is the variable of integration. It is possible to see the same thing if you Author INT(x,@,0,1)

5. Dividing a problem into small independent pieces is usually better. Getting the connection with Gnuplot solved will make it easier once you have the function you want it to render.

8. Try something other than TAYLOR for approximating a DE. Try using those on a simplified problem like dy/dx = -1/x^2 which you can use DSOLVE1(- 1/x^2, -1, x, y, 1, 1) but you can also use TAYLOR and the other tools. See if you can find a tool that will give you a good looking solution on something like that which does not have a nice polynomial solution. I hesitate to give you the answer. Explore what is available and learn how to control what you have.

P.S. Individuals usually think what they know and use is better than what others use. That applies to computer algebra systems, computer languages, brands of cars and lots of other things in the world that I will not mention here. Many people know a little about one item and not so much or even very little about others. You commonly see "I think X is good and have never used or know anything about anything else." Choose your tools and spend lots of time becoming very skilled in using them. That will repay you. If almost everyone around you is using Y then choosing Y may help you share knowledge. But knowing more than one tool and knowing how to learn new ones will also benefit you.

 

[Davide82]

I finally managed to plot the graph!
Here it is in attachment (first plot).
I just used the Runge-Kutta numerical method, i.e.
RK([0.255/x^5·(12·x + 6·x^2 + x^3)·(EXP(-x) - y·(1 + EXP(-x)))], [x, y], [0.5, 0.5], 0.01, 1000)

I only have two questions still:
1) I used x0=0.5 and y0=0.5 but, actually, the initial conditions for the differential equation are y-->0.5 when x-->0. I actually don't know how to set them up properly in the plot: when I make x0 smaller and smaller (let's say 0.1) it plots strange things in the first part of the graph (it is attached as second plot). Have you perhaps some ideas about the reason? Maybe it is caused by errors in Runge-Kutta approximantion? Any ideas how to circumvent them?

2) Since the function I am looking for has an asymptote, I would like to calculate it. Is there a nice way to find it instead of looking at the graph for big values of x-axis?

 

[Bill Simpson]

You say you want an initial condition at x=0, but you have x^5 in your denominator. What do you want .255/0^5 power to do?

For your asymptote, you might consider SUB which subscripts a vector. Could you explore how to use that to extract a single y value from what RK gives you?

 

[Davide82]

I know I have x^5 in the denominator; I only said that when x approaches toward 0 y should be 0.5. Given the fact that I was using a step of 0.01 I believed I could move toward zero at least with x = 0.1 (which is an order of magnitude larger than the step); but in fact it doesn't work. What do you suggest? My aim is only to get a nice plot starting from x=0.1 if possible.

Ok, I try the SUB command. I was wondering if there was an analytical method to find the asymptote of a differential equation.

Thank you.

 

[Bill Simpson]

I don't know what it is about the form of your differential equation, but as I change x from .5 to .4 to .3 you can verify that it begins to blow up and the y values near that smallest x value become extremely large. I have not taken the time to track down exactly why that is happening.

If you can find an analytic solution to your DE then you can use LIM to take a limit and find the asymtote. There are also other more complicated tricks to do that. But when RK only returns you a vector of values and you cannot get an analytic solution I am not aware of any other method of determining the limit.

 

[Davide82]

Thank you Bill.
I am almost satisfied with my plot.
If you have some time and wish to investigate some more it would be nice, but you already helped me a lot and I am very grateful!