3
Attempt to solve Riccati ode y' = f2(x)*y^2+f1(x)*y+f0(x)
7
D Zwillinger, Handbook of Differential Equations, 3rd ed
8
Academic Press, (1997), pp 354-355
10
G M Murphy, Ordinary Differential Equations and Their
11
Solutions, Van Nostrand, 1960, pp 15-23
14
Copyright (C) 2004 David Billinghurst
16
This program is free software; you can redistribute it and/or modify
17
it under the terms of the GNU General Public License as published by
18
the Free Software Foundation; either version 2 of the License, or
19
(at your option) any later version.
21
This program is distributed in the hope that it will be useful,
22
but WITHOUT ANY WARRANTY; without even the implied warranty of
23
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24
GNU General Public License for more details.
26
You should have received a copy of the GNU General Public License
27
along with this program; if not, write to the Free Software
28
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
31
declare(method,special);
32
put('ode1_riccati,001,'version)$
34
ode1_riccati(eq,y,x) := block(
36
de:expand(lhs(eq)-rhs(eq)),
37
%a:coeff(de,'diff(y,x),1),
38
if %a=0 then return(false),
40
f2:-ratsimp(coeff(de,y,2)),
41
if not(freeof(y,f2)) then return(false),
42
if f2=0 then return(false),
43
f1:-expand(ratsimp(coeff(de,y,1))),
44
if not(freeof(y,f1)) then return(false),
45
f0:-expand(ratsimp(de-'diff(y,x)+f2*y^2+f1*y)),
46
if not(freeof(y,f0)) then
48
if not(is(ratsimp(expand(de-'diff(y,x)+f2*y^2+f1*y+f0))=0)) then
51
ode_disp(" is Ricatti equation"),
52
ode_disp2(" f0: ",f0),
53
ode_disp2(" f1: ",f1),
54
ode_disp2(" f2: ",f2),
56
/* Following Murphy, (3-1, p15) see if the equation has the form
57
of the original equation studied by Riccati
67
if ( f1#0 ) then return(false),
68
if (not freeof(x,f2) ) then return(false),
72
if ( ratsimp(f0-c*x^m)#0 ) then return(false),
73
ode_disp(" equation is original Riccati equation"),
74
ode1_riccati_original(b,c,m,y,x)
76
if ( ans#false ) then (
81
/* Perhaps it has the special form Murphy (3-3, p21-22)
83
x*y' -a*y + b*y^2 = c*x^n
85
=> a: f1(x)*x constant
91
if ( not freeof(x,a:ratsimp( f1*x)) ) then return(false),
92
if ( not freeof(x,b:ratsimp(-f2*x)) ) then return(false),
95
/* May want to check c and m */
96
if ( is(ratsimp(f0-c*x^m)#0) ) then return(false),
98
ode_disp(" equation is special Riccati equation"),
99
ode1_riccati_special(a,b,c,n,y,x)
101
if ( ans#false ) then (
106
/* The equation doesn't have a special form, so it is a general
109
ans:ode1_riccati_general(f0,f1,f2,y,x),
110
if ( ans#false ) then (
115
/* Default return value */
119
/* Solve the original Riccati equation y' + b*y^2 = c*x^m
122
ode1_riccati_original(b,c,m,y,x) := block(
124
ode_disp(" -> In ode1_riccati_original"),
129
/* Solve m*(2*k+1)+4*k=0 => k= m/(2*m+4)
131
If k is an integer then the equation is integrable
134
The solution is then found using the transformation y = w/x, giving
136
x*w'(x)-w+b*w^2=c*x^(m+2)
138
which is the special Riccati equation with a=1 and n=m+2
140
if ( not(equal(m,-2)) and integerp(k:m/(2*m+4)) ) then (
141
ode_disp(" Equation is integrable in finite terms"),
142
ode_disp(" Transforming using y=w/x and calling ode1_riccati_special"),
143
ans:ode1_riccati_special(1,b,c,m+2,w,x),
144
ode_disp2(" Solution to transformed equation is ",ans),
145
return(y=ratsimp(rhs(ans)/x))
148
ode1_riccati_original_not_integrable(b,c,m,y,x)
152
/* Solve the original Riccati equation y' + b*y^2 = c*x^m
153
for cases not integrable in finite terms
155
ode1_riccati_original_not_integrable(b,c,m,y,x) := block(
156
if is(equal(m,-2)) then
157
ode1_riccati_original_2(b,c,y,x)
159
ode1_riccati_original_3(b,c,m,y,x)
162
/* Solve the original Riccati equation y' + b*y^2 = c*x^-2
163
- Transform to second order linear ode, Murphy (3-2c, p20-21)
164
- Solve using Murphy A3-250
166
ode1_riccati_original_2(b,c,y,x) := block(
168
ode_disp(" -> In ode1_riccati_original_2"),
169
ode_disp(" Original Riccati equation with m=-2"),
174
/* Let b*y*u(x)=u'(x) so that ode becomes
176
u''(x) - b*c*x^-2*u(x) = 0
178
NOTE: Murphy (p21) has sign of second term wrong.
181
ode_disp2(" r2: ",r2),
184
/* Murphy A3-250-i */
185
ode_disp(" Case i: r2>0"),
188
u:sqrt(x)*(cos(r*log(x))+%c*sin(r*log(x)))
190
else if (r2<0) then (
191
/* Murphy A3-250-ii */
192
ode_disp(" Case ii: r2<0"),
193
r:sqrt(-r2), /* typo in Murphy: missing "-" */
195
u:sqrt(x)*(x^r+%c*x^-r)
197
else if is(equal(r2,0)) then (
198
/* Murphy A3-250-iii */
199
ode_disp(" Case iii: r2=0"),
200
u:sqrt(x)*(1+%c*log(x))
203
error("ode1_riccati_original_2: Impossible case")
206
return(y=ratsimp(diff(u,x)/(b*u)))
209
/* Solve the original Riccati equation y' + b*y^2 = c*x^m for m#-2
210
- Transform to second order linear ode, Murphy (3-2c, p20-21)
211
- Solve using Murphy A3-41
213
ode1_riccati_original_3(b,c,m,y,x) := block(
214
[p:m+2,n:1/(m+2),%c,u,signb,signc],
215
ode_disp(" -> In ode1_riccati_original_3"),
216
ode_disp(" Original Riccati equation with m#-2"),
223
/* Let b*y*u(x)=u'(x) so that ode becomes
225
u''(x) - b*c*x^m*u(x) = 0
227
NOTE: Murphy (p21) has sign of second term wrong.
229
Solution expressed in terms of Bessel functions of order n.
235
if ( (signb='pos and signc='neg) or (signb='neg and signc='pos) ) then
236
if integerp(n) then (
237
u: sqrt(x)*(bessel_j(n,2*sqrt(-b*c)*x^(p/2)/p)
238
+ %c*bessel_y(n,2*sqrt(-b*c)*x^(p/2)/p) )
240
u: sqrt(x)*(bessel_j(n,2*sqrt(-b*c)*x^(p/2)/p)
241
+ %c*bessel_j(-n,2*sqrt(-b*c)*x^(p/2)/p) )
243
else if ((signb='pos and signc='pos) or (signb='neg and signc='neg)) then
244
if integerp(n) then (
245
u: sqrt(x)*(bessel_i(n,2*sqrt(b*c)*x^(p/2)/p)
246
+ %c*bessel_k(n,2*sqrt(b*c)*x^(p/2)/p) )
248
u: sqrt(x)*(bessel_i(n,2*sqrt(b*c)*x^(p/2)/p)
249
+ %c*bessel_k(-n,2*sqrt(b*c)*x^(p/2)/p) )
252
/* b and c are non-zero constants, so this is an error */
253
error("ode_riccati_original_3: Impossible case has just happened"),
254
return(y=ratsimp(diff(u,x)/(b*u)))
257
/* Solve the special Riccati equation x*y' -a*y + b*y^2 = c*x^n
260
ode1_riccati_special(a,b,c,n,y,x) := block(
261
[k,s,u,%c,signb,signc],
262
ode_disp(" -> In ode1_riccati_special"),
268
/* Certain cases are integrable. */
271
Equation can be made exact using integrating factor x^(a-1)
274
if ( is(equal(n,2*a)) ) then (
275
ode_disp(" Case (a.i)"),
276
return(ode1_riccati_special_i(a,b,c,n,y,x))
279
/* Case (a.ii) (n-2*a)/(2*n) a positive integer */
280
else if ( n#0 and integerp(k:(n-2*a)/(2*n)) and k>0 ) then (
281
ode_disp2(" Case (a.ii) with k = ",k),
283
s:rhs(ode1_riccati_special_i(n/2,c,b,n,y,x))
285
s:rhs(ode1_riccati_special_i(n/2,b,c,n,y,x)),
286
for i:(k-1) thru 1 step -1 do (
295
/* Case (a.iii) (n+2*a)/(2*n) a positive integer */
296
else if ( n#0 and integerp(k:(n+2*a)/(2*n)) and k>0 ) then (
297
ode_disp2(" Case (a.iii) with k = ",k),
299
s:rhs(ode1_riccati_special_i(n/2,c,b,n,y,x))
301
s:rhs(ode1_riccati_special_i(n/2,b,c,n,y,x)),
302
for i:(k-1) thru 1 step -1 do (
311
/* Not integrable in finite terms. For a=0 we have
314
Let y = x*u'(x)/(b*u(x)) => x^2*u'' + x*u' - b*c*x^n*u = 0
316
This is Murphy A3.202 or Abramowitz and Stegun 9.1.53
318
else if is(equal(a,0)) then (
319
/* Can never have n=0 so division is always safe */
320
ode_disp(" Case (c.i)"),
325
if ( (signb='pos and signc='neg) or (signb='neg and signc='pos) ) then
326
u: bessel_j(0,2*sqrt(-b*c)*x^(n/2)/n)
327
+ %c*bessel_y(0,2*sqrt(-b*c)*x^(n/2)/n)
329
else if ((signb='pos and signc='pos) or (signb='neg and signc='neg)) then
330
u: bessel_i(0,2*sqrt(b*c)*x^(n/2)/n)
331
+ %c*bessel_k(0,2*sqrt(b*c)*x^(n/2)/n)
333
/* b and c are non-zero constants, so this is an error */
334
error("ode_riccati_special: Impossible case has just happened"),
336
return(y=ratsimp(x*diff(u,x)/(b*u)))
338
/* For a#0, transform using y=z*u(z) with z=x^a
339
to give u' + (b/a)*u^2 = (c/a)*z^((n-2*a)/a)
340
which is the original Riccati equation
343
ode_disp(" Case (c.ii)"),
344
s:ode1_riccati_original_not_integrable(b/a,c/a,(n-2*a)/a,u,z),
345
ode_disp2(" Solution of transformed eqn is ",s),
346
return(y=ratsimp(subst(x^a,z,z*rhs(s))))
350
/* Solve the special Riccati equation x*y' -a*y + b*y^2 = c*x^n
351
for the case n=2*a. Murphy (3-3, p21-22).
353
Note: Signs changed from Murphy in cases a.i.1 and a.i.3.
355
ode1_riccati_special_i(a,b,c,n,y,x) := block(
357
ode_disp(" -> In ode1_riccati_special_i"),
363
if not(equal(n,2*a)) then error("ode1_riccati_special_i: n#2*a"),
366
Equation can be made exact using integrating factor x^(a-1)
367
and integrated to give solution(s) below.
372
if ( signb='pos and signc='pos ) then (
373
/* Murphy has the sign of the solution wrong */
374
ode_disp(" Case (a.i.1) b*c>0, b>0 and c>0"),
375
return(y=sqrt(c/b)*x^a*tanh(sqrt(b*c)*x^a/a+%c))
377
else if ( signb='neg and signc='neg ) then (
378
ode_disp(" Case (a.i.2) b*c>0, b<0 and c<0"),
379
return(y=sqrt(c/b)*x^a*tanh(%c-sqrt(b*c)*x^a/a))
382
else if ( signb='pos and signc='neg ) then (
383
ode_disp(" Case (a.i.3) b*c<0, b>0 and c<0"),
384
return(y=sqrt(-c/b)*x^a*tan(%c-sqrt(-b*c)*x^a/a))
386
else if ( signb='neg and signc='pos ) then (
387
ode_disp(" Case (a.i.4) b*c<0, b<0 and c>0"),
388
/* Murphy has the sign of the solution wrong */
389
return(y=sqrt(-c/b)*x^a*tan(sqrt(-b*c)*x^a/a+%c))
392
/* b and c are non-zero constants, so this is an error */
393
error("ode_riccati_special_i: Impossible case has just happened")
396
/* The equation doesn't have a special form, so it is a generalized
397
Riccati equation. Try transforming it to a linear second order
400
Substitute y = -z'/(z*f2)
402
=> f2*z''-(f2'+f1*f2)z'+f2^2*f0*z=0
404
Solve this second order linear ode for z. The solution has form
405
z=%k1*f+%k2*g, with two constants, but a first order ode only has
406
one constant %c. Without loss of generality take %k1=1 and %k2=%c
408
y = -z'/(z*f2) = -(f'+%c*g')/((f+%c*g)*f2)
410
ode1_riccati_general(f0,f1,f2,y,x) := block(
411
[de,z,ans,%c,%k1,%k2],
412
ode_disp(" Transforming to 2nd order ode"),
413
de: f2*'diff(z,x,2)-(diff(f2,x)+f1*f2)*'diff(z,x)+f2^2*f0*z=0,
414
if get('contrib_ode,'verbose) then disp(de),
416
ans:contrib_ode(de,z,x),
417
ode_disp(" with solution"),
418
if get('contrib_ode,'verbose) then disp(ans),
420
if is(ans=false) then (
421
ode_disp(" Cannot solve 2nd order ode"),
424
else if (lhs(ans[1])#z) then (
425
/* give up on parametric solutions for z */
426
ode_disp(" 2nd order ode has parametric solution"),
427
ode_disp(" giving up"),
431
ode_disp(" solution found"),
433
ans:rhs(ans[1]), ans:subst(1,%K1,ans), ans:subst(%C,%K2,ans),
434
return(y=ratsimp(-diff(ans,x)/(ans*f2)))
added
removed
share/contrib/diffequations/ode1_riccati.mac
