deSolve how to do this on Prime



#18

desolve example: how do I do this on Prime?

I can't figure out how to execute desolve on an example from the NSpire desolve entry in their refernce guide. The example is: deSolve(y''-2*y'/x+(9+2/x^2)*y=x*e^x and y(pi/6)=0 and y(pi/3)=0,x,y)

I used the y'' example as a template from the onboard help. I get empty brackets in CAS and syntax error in home.


#19

I am trying to find out, but I can't solve this in any system I've tried yet including mathematica. Anyone with more experience with ode solving want to chime in with is needed there?

Result spit out by the nspire (shown in 2d at wolfram) is: wolfram 2d display

TW

Edited: 21 Oct 2013, 3:38 p.m.


#20

Well, even the good old TI-92+ has no problems to solve this ODE:

y=C1*x*cos(3x)+C2*x*sin(3x)+x*e^x/10 (general solution)

And with the 2 initial conditions: C1=e^(pi/3)/10, C2=-e^(pi/6)/10

It seems that XCAS could learn something from the old Derive CAS ... ;-)

Franz


#21

And mathematica apparently? Don't think so.

Question still stands as to what is happening there since no other tool I've tried likes that input or other variants I've tried and gives *any* result - let alone a similar one.

TW

Edited: 21 Oct 2013, 4:50 p.m.


#22

Quote:
Question still stands as to what is happening there since no other tool I've tried likes that input or other variants I've tried and gives *any* result - let alone a similar one.

I've now checked the solution and it's definitely correct.

So the question is not "why Derive or TI-92 solve it" but "why do the other CAS not solve it"?

Franz


#23

I believe that was the original question.

TW

Edited: 21 Oct 2013, 5:53 p.m.

#24

Quote:
It seems that XCAS could learn something from the old Derive CAS ... ;-)

You will always find some very specific problem where a miraculous conjunction of parameters value leads to a simple answer. You could even do that by choosing the solution and generating coefficients so that it solves the corresponding equation, then you put it in the examples of the manual, and then you tell to the world look I'm solving this equation, other don't, I'm better.
I'm not saying derive does that but the fact is that generic linear differential equations with non constant coefficients do not have solutions expressed in terms of elementary functions, that's why we define Bessel functions for example. I never bothered to check for very specific conjunction of parameters value where they have. Now of course if someone knows a nice simple algorithm to do that, I'd be happy to implement it!

#25

Quote:
Now of course if someone knows a nice simple algorithm to do that, I'd be happy to implement it!

Well, no problem! ;-)

The method for linear 2nd order ODEs with constant coefficients can easily

be extended for 2 other types of (similar) ODEs - one an even more general

form and the other one a special type (Euler DE).

I've tried these 2 other types in XCAS and got an "Unable to solve ..." message,

so here's a 'simple algorithm' for implementing in XCAS and HP-Prime.

It's the method I've used in my own ODE-Solver package for the TI-92, written

long time ago because the first TI-92 versions had no ODE solving commands,

this has been implemented later in the TI-92+.

(in pseudo-code, because not everyone knows the TI-92 language)

y''=u*y'+v*y+w  (with u,v,w functions of x)

Type 1 (constant coefficients): u=const and v=const

Type 2 (more general form of type 1):
u and v not constant, but u^2/4-u'/2+v=const

Type 3 (special form of type 1, Euler differential equation):
u=a/x and v=b/x^2 (with a,b=const)

Assuming you have a linear 2nd order ODE of the form y''=u*y'+v*y+w,
and the coefficient terms u, v and w are already assigned, here's the
algorithm for the 3 types mentioned above:

k:=u^2/4-u'/2+v
if k=const or k*x^2=const then
if k=const then
s:=x
t:=e^(int(u,x)/2)
else
u:=u*x+1
k:=u^2/4+v*x^2
s:=ln(x)
t:=x^(u/2)
endif
if k=0 then
u:=t*s
v:=t
elseif k>0 then
u:=t*e^(sqrt(k)*s)
v:=t*e^(-sqrt(k)*s)
else
u:=t*cos(sqrt(-k)*s)
v:=t*sin(sqrt(-k)*s)
endif
w:=w/(u*v'-v*u')
w:=v*int(u*w,x)-u*int(v*w,x)
solution: y=c1*u+c2*v+w
endif

Franz


#26

Great! I'll do it, many thanks!


#27

I've now tried a few other ODEs in XCAS, and found some cases where XCAS crashes (i.e. simply exits) - here's one example:

deSolve(y''=x*y'+x^2*(y')^2)

I wanted to test if XCAS can solve Liouville's ODEs of the form
y''=u(x)*y'+v(y)*y'^2, but I simply overlooked that the coefficient v(y) of y'^2 should be a function of y (instead of x).

So my example above was of course no Liouville ODE but even a simpler one (called 'autonom' in German, which means y''=f(x,y') or y''=f(y,y'), i.e. either y or x does not exist in the RHS) which is usually solved by reducing it to a 1st order ODE.

Maybe my example above is not solvable by XCAS, but of course it should not crash!

Franz

Edited: 22 Oct 2013, 8:32 a.m.


#28

Just tried with the latest source, I can't reproduce a crash. I get
[integrate((c_0+integrate(-x^2*exp(x^2/2),x))^-1*exp(x^2/2),x)+c_1]
By the way, in France we call autonom an equation independant of the time x, for this one it's not the case, we call it "incomplète en y" (don't know how this is translated in English).


#29

Quote:
Just tried with the latest source, I can't reproduce a crash. I get
[integrate((c_0+integrate(-x^2*exp(x^2/2),x))^-1*exp(x^2/2),x)+c_1]

I've also used the latest XCAS Windows version from yesterday (xcas 1.1.0-20), and and xcas just exits with this example.

It shows a few lines of error messages, but so quickly that I can't read them.

Franz


#30

Very strange, as I said it works on my devel box (32 bits linux) but I could reproduce the segfault under linux 64 bits, the messages are
Temporary replacing surd/NTHROOT by fractional powers
Seems to be an infinite loop calling integration, will investigate tomorrow.
Edit:
Well it seems gcc does not always call the same routine when -O2 optimization is on or not, which explains the different behavior. Changing one of the homonyms seems to solve the problem, I'll need to make a 1.1.0-21 update for xcas.


Edited: 22 Oct 2013, 3:16 p.m.


#31

Quote:
I'll need to make a 1.1.0-21 update for xcas.

Yep, yesterdays Xcas update fixed this problem and now it can also solve the 2 additional ODE types for which I've posted the solution method above - many thanks for this! :-)

But I'm still a bit disappointed about Xcas' solving capabilities of 2nd order ODEs.

I've tested different kinds of (2nd order) ODEs and it seems that Xcas can only solve 3 types:

1) linear with constant coefficient: y''=a*y'+b*y+s(x) (a,b=const)

2) linear of the form we discussed above (since the last Xcas version)

3) independent of y: y''=f(x,y')

Especially for the 3rd type I wonder why you didn't also implement the type "independent of x", i.e. ODEs of the form y''=f(y,y') ?

This "free of x" type is also very easy to solve and the method is quite similar to the "free of y" type - the only diffenence is that in the substitution u=y' (which requires solving 2 first-order ODEs) here u is a function of y (instead of x), and thus u' has to be built differently (u'=du/dy*dy/dx=u'(y)*y'(x)=u'(y)*u).

Here's the usual way such ODEs y''=f(y,y') are solved:

Given an ODE of the form y''=f(y,y'):

1) substitute y'->u, y''->u'*u
2) leads to u'*u=f(y,u) or u'=f(y,u)/u
3) deSolve(u'=f(y,u)/u,y,u) gives intermediate solution u=u(y)
4) backsubstitution leads to y'=u(y)
5) deSolve(y'=u(y),x,y) gives final solution y=y(x).

Implementing this method would be a big enhancement, because many ODEs are exactly of this type (BTW, also TI-92+ and TI-Nspire can solve this type!).

In my own ODE package for the old TI-92 I even have 2 further 2nd order ODEs implemented (the Liouville-DE y''=u(x)*y'+v(y)*y'^2 and 'exact' 2nd order ODEs), but these 2 are rather seldom used (and a bit more complicated to solve), so I think you won't be interested in them.

But this "free of x" type is indeed very important IMO ...

Franz

Edited: 24 Oct 2013, 4:19 p.m.

#32

So out of curiosity, what happens when you replace any of those constants with *any* slightly different numbers, or a variable on the nspire?

For example, deSolve(y''-((a*y')/(x))+(b+((c)/(x^(2))))*y=x*e^(x) and y(((pi)/(6)))=0 and y(((pi)/(3)))=0,x,y)

TW

Edited: 22 Oct 2013, 12:18 a.m.


#33

The NSpire comes up with y''=(x^3 * e^x - (b*x^2)*y+ a*x*y' - c*y)/x^2.

Substituting arbitrary numbers into a,b,c in the original example mostly results in an answer in the form of y''=

Changing initial conditions is ok...comes up with an answer for y


#34

Quote:
Substituting arbitrary numbers into a,b,c in the original example mostly results in an answer in the form of y''=

It is only solvable if a^2+2a-4c=0 (or if b=0), because only then it's the special type 2 (or type 3) that I've mentioned in my other posting above.

(b is arbitrary as long as it is constant!)

Franz


Edited: 22 Oct 2013, 10:33 a.m.


Possibly Related Threads...
Thread Author Replies Views Last Post
  desolve example Richard Berler 2 369 10-21-2013, 11:09 AM
Last Post: Richard Berler
  desolve Richard Berler 1 274 10-02-2013, 01:59 PM
Last Post: parisse
  HP50g DESOLVE Help Chuck 9 688 04-26-2007, 07:55 PM
Last Post: Chuck

Forum Jump: