Post Reply 
deTaylor
05-22-2015, 01:24 PM (This post was last modified: 05-22-2015 01:39 PM by fhub.)
Post: #1
deTaylor
Here's a small CAS function deTaylor for the HP-Prime (similiar to deSolve) which I've originally written for Xcas.
It approximately solves 1st- and 2nd-order differential equations (with initial conditions) as n-th degree Taylor polynomial.

1st order y'=f(x,y) with y(x0)=y0:
deTaylor(f,[x,y],[x0,y0],n)

2nd order y''=f(x,y,y') with y(x0)=y0 and y'(x0)=y0':
deTaylor(f,[x,y,z],[x0,y0,z0],n)
(I'm using z instead of y', because y' can't be used as input)

Here's the function definition:
(you can directly copy&paste it into the Prime-emulator commandline in CAS-mode)
Code:

deTaylor(f,v,c,n):=
BEGIN
  local k,yk,sol;
  yk:=v(2); sol:=c(2);
  IF dim(v)<2 or dim(v)>3 or dim(v)<>dim(c) THEN
    RETURN("Argument error!");
  END;
  FOR k FROM 1 TO n STEP 1 DO
    IF dim(v)=2 THEN
      yk:=diff(yk,v(1))+f*diff(yk,v(2));
      sol+=subst(yk,[v(1)=c(1),v(2)=c(2)])/k!*(v(1)-c(1))^k;
    ELSE
      yk:=diff(yk,v(1))+v(3)*diff(yk,v(2))+f*diff(yk,v(3));
      sol+=subst(yk,[v(1)=c(1),v(2)=c(2),v(3)=c(3)])/k!*(v(1)-c(1))^k;
    END;
  END;
  expand(simplify(sol));
END

Example 1: y'=x*y^2+1 with y(0)=1 (6th-degree approximation):
deTaylor(x*y^2+1,[x,y],[0,1],6)

Example 2: y''=x*y*y' with y(1)=2 and y'(1)=3 (5th degree):
deTaylor(x*y*z,[x,y,z],[1,2,3],5)

Maybe it's useful for someone, Wink
Franz
Visit this user's website Find all posts by this user
Quote this message in a reply
05-22-2015, 06:10 PM (This post was last modified: 05-22-2015 06:10 PM by salvomic.)
Post: #2
RE: deTaylor
thank you!

please help to control with an differential equation still not solvable in Prime (but now ok in XCas):
y'=(x+y)^2
Is it correct to input
deTaylor((x+y)^2, [x,y], [0,0], 7) ?

I get (17/315)x^7+(2/15)x^5+(⅓)x^3
The general solution of equation (XCas) is TAN(x-c)-x
With Taylor I've taylor(TAN(x)-x), x, 6) = (⅓)x^3+(2/15)x^5+x^7+o(x)
It should be ok, isn't it? Smile

Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
05-22-2015, 08:27 PM
Post: #3
RE: deTaylor
(05-22-2015 06:10 PM)salvomic Wrote:  thank you!

please help to control with an differential equation still not solvable in Prime (but now ok in XCas):
y'=(x+y)^2
Is it correct to input
deTaylor((x+y)^2, [x,y], [0,0], 7) ?

I get (17/315)x^7+(2/15)x^5+(⅓)x^3
The general solution of equation (XCas) is TAN(x-c)-x
With Taylor I've taylor(TAN(x)-x), x, 6) = (⅓)x^3+(2/15)x^5+x^7+o(x)
It should be ok, isn't it? Smile
Not exactly, you should enter order 7 (not 6), then you get the same result as with deTaylor (despite of the error term "x^8*order_size(x)"):

taylor(tan(x)-x), x=0, 7)

Franz
Visit this user's website Find all posts by this user
Quote this message in a reply
05-22-2015, 08:31 PM
Post: #4
RE: deTaylor
(05-22-2015 08:27 PM)fhub Wrote:  Not exactly, you should enter order 7 (not 6), then you get the same result as with deTaylor (despite of the error term "x^8*order_size(x)"):
taylor(tan(x)-x), x=0, 7)
Franz

right!
thank you.
Parisse has already added the solution for y'=(x+y)^2 in XCas, but it was after the last FW was ready. I hope this will be added in the next firmware.
In the meantime your deTaylor is helpful also for this equation ;-)

Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
05-22-2015, 09:44 PM (This post was last modified: 05-22-2015 09:48 PM by fhub.)
Post: #5
RE: deTaylor
(05-22-2015 08:31 PM)salvomic Wrote:  Parisse has already added the solution for y'=(x+y)^2 in XCas, but it was after the last FW was ready. I hope this will be added in the next firmware.
In the meantime your deTaylor is helpful also for this equation ;-)
Well, in the meantime you could also use the following function: Smile
Code:

deLinSubst(f,x,y):=
BEGIN
  local k;
  k:=simplify(diff(f,x)/diff(f,y));
  IF diff(k,x)<>0 or diff(k,y)<>0 THEN
    RETURN("No linear substitution!");
  END;
  RETURN(solve(subst(integrate(1/(subst(f,y=y-k*x)+k),y),y=y+k*x)=x+G_0,y));
END

Works only for ODEs of the type y'=f(a*x+b*y+c) with a,b,c=const.
For your example y'=(x+y)^2 just enter:
deLinSubst((x+y)^2,x,y)
(enter only the RHS, not y'=...)

In this simple form you get a solution only if Xcas can 'solve' the equation for y, but it could easily be modified to return at least an implicit solution if the equation is unsolvable.

Franz
Visit this user's website Find all posts by this user
Quote this message in a reply
05-22-2015, 10:13 PM
Post: #6
RE: deTaylor
(05-22-2015 09:44 PM)fhub Wrote:  Well, in the meantime you could also use the following function: Smile
Code:

deLinSubst(f,x,y):=
...

Works only for ODEs of the type y'=f(a*x+b*y+c) with a,b,c=const.
For your example y'=(x+y)^2 just enter:
deLinSubst((x+y)^2,x,y)
(enter only the RHS, not y'=...)

In this simple form you get a solution only if Xcas can 'solve' the equation for y, but it could easily be modified to return at least an implicit solution if the equation is unsolvable.

Franz

well done, Franz!
it works: gives the true (general) solution of ODE y'=(x+y)^2 -> TAN(G_0+x)-x

Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 3 Guest(s)