#192 Modified Bessel differential equation

[Evert Provoost]

Dear all,

I was trying to see whether Maxima could serve as an alternative to our current tool used in (engineering) calculus courses, partially to settle some friendly banter with my supervisor about FOSS software and partially because the license fee for the currently used package is becoming fairly steep . As I expected from my experience through the years, Maxima does very well. Especially cases where it is more 'pedantic' are didactically valuable as our current tool often gives a result which is incorrect for some edge case.

One introductory example where it does give some less than ideal results is the modified Bessel differential equation:

mbdv: x^2*'diff(y, x, 2) + x*'diff(y, x) - (x^2 + n^2)*y = 0$
ode2(mbdv, y, x);

which gives false as output. contrib_ode does give a result, but it appears somewhat convoluted. Would it be possible to add the textbook solution using the modified Bessel functions?

Already much thanks for making this wonderful piece of software ever the more powerful! (If I were to come across further hiccups in this experiment I'll let you know.)

[Evert Provoost]

I seem to have clicked the wrong button, this is of course a Feature Request, not a Bug... Can this be moved, or do I reopen this under the right category?

[Robert Dodier]

Ticket moved from /p/maxima/bugs/4389/

[Robert Dodier]

Evert, I've moved this ticket to the feature requests. However, since some discussion seems appropriate, I'd like to suggest that you bring up this topic on the maxima-discuss mailing list. If you are not already subscribed, you can subscribe on this web page: https://sourceforge.net/projects/maxima/lists/maxima-discuss

[David Billinghurst]

The function of interest is maxima code - the nine-line function bessel2 in file share/diffequations/ode2.mac.

[David Billinghurst]

This looks "trivial" to implement: clone the function bessel2, flip a couple of signs, change the Bessel function names, call the new function and fix the testsuite. I may have time to investigate later this week.

[David Billinghurst]

Here is an initial implementation. Works for simple cases, but not for n=3/2. Enough for tonight.

Here is a patch to ode2.mac. Introduce new function modifiedbessel2.mac and call it.

Evert: I can send you a complete patched ode2.mac if you aren't familiar with the patch utility.

$ diff -c /usr/local/src/maxima/share/diffequations/ode2.mac ode2.mac
*** /usr/local/src/maxima/share/diffequations/ode2.mac  Mon Sep 12 21:21:58 2016
--- ode2.mac    Tue Oct 22 21:54:59 2024
***************
*** 11,17 ****
  declare_translated(boundtest,noteqn,nlxy,nly,nlx,xcc2,bessel2,euler2,
     pttest,exact2,cc2,genhom,solvebernoulli,solvehom,integfactor,exact,
     separable,solvelnr,solve1,linear2,reduce,hom2,pr2,varp,desimp,failure,
!    ode1a,ftest,ode2a));

  block(local(eq,yold,x),
  ode2(eq,yold,x):=block([derivsubst:false,ynew],
--- 11,17 ----
  declare_translated(boundtest,noteqn,nlxy,nly,nlx,xcc2,bessel2,euler2,
     pttest,exact2,cc2,genhom,solvebernoulli,solvehom,integfactor,exact,
     separable,solvelnr,solve1,linear2,reduce,hom2,pr2,varp,desimp,failure,
!    ode1a,ftest,ode2a,modifiedbessel2));

  block(local(eq,yold,x),
  ode2(eq,yold,x):=block([derivsubst:false,ynew],
***************
*** 172,179 ****
     if ftest(cc2(ap,aq,y,x)) then return(%q%),
     if ftest(exact2(a1,a2,a3)) then return(%q%),
     if (pt: pttest(ap)) = false then go(end),
!    if ftest(euler2(ap,aq)) then return(%q%),
!    if ftest(bessel2(ap,aq)) then return(%q%),
   end,
     if ftest(xcc2(ap,aq)) then return(%q%) else return(false)))$

--- 172,180 ----
     if ftest(cc2(ap,aq,y,x)) then return(%q%),
     if ftest(exact2(a1,a2,a3)) then return(%q%),
     if (pt: pttest(ap)) = false then go(end),
!    if ftest(euler2(ap,aq,pt)) then return(%q%),
!    if ftest(bessel2(ap,aq,pt)) then return(%q%),
!    if ftest(modifiedbessel2(ap,aq,pt)) then return(%q%),
   end,
     if ftest(xcc2(ap,aq)) then return(%q%) else return(false)))$

***************
*** 293,300 ****
     a3: coeff(a1,x,0),
     if not(a1 = a2*x + a3) then return(false) else return(-a3/a2)))$

! block(local(a,b),
! euler2(a,b):=block([dc,rp,ip,alpha,beta,sign,radexpand:false,%k1,%k2,pos,zero],
     if not(freeof(x,y,beta: ratsimp(b*(x-pt)^2))) then return(false),
     method: 'euler, alpha: a*(x-pt),
     dc: ratsimp((alpha-1)^2 - 4*beta),
--- 294,301 ----
     a3: coeff(a1,x,0),
     if not(a1 = a2*x + a3) then return(false) else return(-a3/a2)))$

! block(local(a,b,pt),
! euler2(a,b,pt):=block([dc,rp,ip,alpha,beta,sign,radexpand:false,%k1,%k2,pos,zero],
     if not(freeof(x,y,beta: ratsimp(b*(x-pt)^2))) then return(false),
     method: 'euler, alpha: a*(x-pt),
     dc: ratsimp((alpha-1)^2 - 4*beta),
***************
*** 306,313 ****
     dc: -dc, ip: sqrt(dc)/2,
     return(y = (x-pt)^rp * (%k1*sin(ip*log(x-pt)) + %k2*cos(ip*log(x-pt))))))$

! block(local(a,b),
! bessel2(a,b):=block([nu,b1,intp,radexpand:'all,%k1,%k2],
     if not(freeof(x,y,b1: ratsimp((1-b)*(x-pt)^2))) then return(false),
     if ratsimp(a*(x-pt)) # 1 then return(false),
     nu: sqrt(b1), method: 'bessel,
--- 307,314 ----
     dc: -dc, ip: sqrt(dc)/2,
     return(y = (x-pt)^rp * (%k1*sin(ip*log(x-pt)) + %k2*cos(ip*log(x-pt))))))$

! block(local(a,b,pt),
! bessel2(a,b,pt):=block([nu,b1,intp,radexpand:'all,%k1,%k2],
     if not(freeof(x,y,b1: ratsimp((1-b)*(x-pt)^2))) then return(false),
     if ratsimp(a*(x-pt)) # 1 then return(false),
     nu: sqrt(b1), method: 'bessel,
***************
*** 316,321 ****
--- 317,332 ----
     if intp = 'yes then return(y = %k1*bessel_j(nu,x-pt) + %k2*bessel_y(nu,x-pt)),
     return(y = %k1*bessel_j(nu,x-pt) + %k2*bessel_j(-nu,x-pt))))$

+ block(local(a,b,pt),
+ modifiedbessel2(a,b,pt):=block([nu,b1,intp,radexpand:'all,%k1,%k2],
+    if not(freeof(x,y,b1: ratsimp(-(1+b)*(x-pt)^2))) then return(false),
+    if ratsimp(a*(x-pt)) # 1 then return(false),
+    nu: sqrt(b1), method: 'modifiedbessel,
+    if nu = 1/2 then return(y = (%k1*sinh(x-pt) + %k2*exp(-(x-pt)))/sqrt(x-pt)),
+    intp:askinteger(nu),
+    if intp = 'yes then return(y = %k1*bessel_i(nu,x-pt) + %k2*bessel_k(nu,x-pt)),
+    return(y = %k1*bessel_j(nu,x-pt) + %k2*bessel_j(-nu,x-pt))))$
+
  block(local(soln,xc,yc),
  ic1(soln,xc,yc):=
     block([%c],

Some simple tests. All good, except n=3/2

(%i1) batch("modifiedbessel.mac");

read and interpret /home/dabilling/maxima/ode2/ModifiedBessel/modifiedbessel.mac
(%i2) display2d:false
(%i3) load("./ode2.mac")
(%i4) load('contrib_ode)
(%i5) mbde:x^2*'diff(y,x,2)+x*'diff(y,x)+(-(x^2+n^2))*y = 0
(%i6) ev(eq:mbde,n = 1)
(%o6) x^2*'diff(y,x,2)+x*'diff(y,x,1)+(-x^2-1)*y = 0
(%i7) s:ode2(eq,y,x)
(%o7) y = bessel_k(1,x)*%k2+bessel_i(1,x)*%k1
(%i8) [method,r:ode_check(eq,s)]
(%o8) [modifiedbessel,0]
(%i9) ev(eq:mbde,n = 2)
(%o9) x^2*'diff(y,x,2)+x*'diff(y,x,1)+(-x^2-4)*y = 0
(%i10) s:ode2(eq,y,x)
(%o10) y = bessel_k(2,x)*%k2+bessel_i(2,x)*%k1
(%i11) [method,r:ode_check(eq,s)]
(%o11) [modifiedbessel,0]
(%i12) ev(eq:mbde,n = 3)
(%o12) x^2*'diff(y,x,2)+x*'diff(y,x,1)+(-x^2-9)*y = 0
(%i13) s:ode2(eq,y,x)
(%o13) y = bessel_k(3,x)*%k2+bessel_i(3,x)*%k1
(%i14) [method,r:ode_check(eq,s)]
(%o14) [modifiedbessel,0]
(%i15) ev(eq:mbde,n = 1/2)
(%o15) x^2*'diff(y,x,2)+x*'diff(y,x,1)+(-x^2-1/4)*y = 0
(%i16) s:ode2(eq,y,x)
(%o16) y = (%k1*sinh(x)+%k2*%e^-x)/sqrt(x)
(%i17) [method,r:ode_check(eq,s)]
(%o17) [modifiedbessel,0]
(%i18) ev(eq:mbde,n = 3/2)
(%o18) x^2*'diff(y,x,2)+x*'diff(y,x,1)+(-x^2-9/4)*y = 0
(%i19) s:ode2(eq,y,x)
(%o19) y = bessel_j(-(3/2),x)*%k2+bessel_j(3/2,x)*%k1
(%i20) [method,r:ode_check(eq,s)]
(%o20) [modifiedbessel,
        (bessel_j(1/2,x)*%k2*x^2)/4-(3*bessel_j(-(3/2),x)*%k2*x^2)/2
                                   +(bessel_j(-(7/2),x)*%k2*x^2)/4
                                   +(bessel_j(7/2,x)*%k1*x^2)/4
                                   -(3*bessel_j(3/2,x)*%k1*x^2)/2
                                   +(bessel_j(-(1/2),x)*%k1*x^2)/4
                                   -(bessel_j(-(1/2),x)*%k2*x)/2
                                   +(bessel_j(-(5/2),x)*%k2*x)/2
                                   -(bessel_j(5/2,x)*%k1*x)/2
                                   +(bessel_j(1/2,x)*%k1*x)/2
                                   -(9*bessel_j(-(3/2),x)*%k2)/4
                                   -(9*bessel_j(3/2,x)*%k1)/4]
(%i21) declare(n,integer)
(%o21) done
(%i22) s:ode2(mbde,y,x)
(%o22) y = %k2*bessel_k(n,x)+%k1*bessel_i(n,x)
(%i23) [method,r:ode_check(mbde,s)]
(%o23) [modifiedbessel,0]

and the test file using a patched local copy of ode2.mac

/* modifiedbessel.mac */

display2d:false$
load("./ode2.mac")$
load('contrib_ode)$

mbde: x^2*'diff(y, x, 2) + x*'diff(y, x) - (x^2 + n^2)*y = 0$

eq: mbde, n=1;
s:ode2(eq,y,x);
[method,r:ode_check(eq,s)];

eq: mbde, n=2;
s:ode2(eq,y,x);
[method,r:ode_check(eq,s)];

eq: mbde, n=3;
s:ode2(eq,y,x);
[method,r:ode_check(eq,s)];

eq:mbde,n=1/2;
s:ode2(eq,y,x);
[method,r:ode_check(eq,s)];

eq:mbde,n=3/2;
s:ode2(eq,y,x);
[method,r:ode_check(eq,s)];

declare(n,integer);
s:ode2(mbde,y,x);
[method,r:ode_check(mbde,s)];

[David Billinghurst]

Ah. I forgot I could add attachments.

Edit. ode2.mac updated with one-line patch below.

[David Billinghurst]

The last line of function modifiedbessel2 should read

return(y = %k1*bessel_i(nu,x-pt) + %k2*bessel_k(nu,x-pt))))$

Case n=3/2 now works. It isn't simplified to zero by ode_check, but the result simplifies to zero numerically.

(%i32) eq:mbde,n=3/2;
(%o32) x^2*'diff(y,x,2)+x*'diff(y,x,1)+(-x^2-9/4)*y = 0
(%i33) s:ode2(eq,y,x);
(%o33) y = bessel_k(3/2,x)*%k2+bessel_i(3/2,x)*%k1
(%i34) [method,r:ode_check(eq,s)];
(%o34) [modifiedbessel,
        (bessel_k(7/2,x)*%k2*x^2)/4-(bessel_k(3/2,x)*%k2*x^2)/2
                                   +(bessel_k(1/2,x)*%k2*x^2)/4
                                   +(bessel_i(7/2,x)*%k1*x^2)/4
                                   -(bessel_i(3/2,x)*%k1*x^2)/2
                                   +(bessel_i(-(1/2),x)*%k1*x^2)/4
                                   -(bessel_k(5/2,x)*%k2*x)/2
                                   -(bessel_k(1/2,x)*%k2*x)/2
                                   +(bessel_i(5/2,x)*%k1*x)/2
                                   +(bessel_i(1/2,x)*%k1*x)/2
                                   -(9*bessel_k(3/2,x)*%k2)/4
                                   -(9*bessel_i(3/2,x)*%k1)/4]
(%i35) r,x=1.0;
(%o35) 1.8865117801247777e-16*%k1
(%i36) r,x=2.0;
(%o36) 2.1510571102112408e-15*%k1
(%i37) r,x=3.0;
(%o37) -(1.1102230246251565e-16*%k2)-5.728750807065808e-14*%k1

[David Billinghurst]

The last line of function modifiedbessel2 should read

return(y = %k1*bessel_i(nu,x-pt) + %k2*bessel_k(nu,x-pt))))$

Case n=3/2 now works. It isn't simplified to zero by ode_check, but the result simplifies to zero numerically.

This ode_check issue is [bug:#4392]

(%i32) eq:mbde,n=3/2;
(%o32) x^2*'diff(y,x,2)+x*'diff(y,x,1)+(-x^2-9/4)*y = 0
(%i33) s:ode2(eq,y,x);
(%o33) y = bessel_k(3/2,x)*%k2+bessel_i(3/2,x)*%k1
(%i34) [method,r:ode_check(eq,s)];
(%o34) [modifiedbessel,
        (bessel_k(7/2,x)*%k2*x^2)/4-(bessel_k(3/2,x)*%k2*x^2)/2
                                   +(bessel_k(1/2,x)*%k2*x^2)/4
                                   +(bessel_i(7/2,x)*%k1*x^2)/4
                                   -(bessel_i(3/2,x)*%k1*x^2)/2
                                   +(bessel_i(-(1/2),x)*%k1*x^2)/4
                                   -(bessel_k(5/2,x)*%k2*x)/2
                                   -(bessel_k(1/2,x)*%k2*x)/2
                                   +(bessel_i(5/2,x)*%k1*x)/2
                                   +(bessel_i(1/2,x)*%k1*x)/2
                                   -(9*bessel_k(3/2,x)*%k2)/4
                                   -(9*bessel_i(3/2,x)*%k1)/4]
(%i35) r,x=1.0;
(%o35) 1.8865117801247777e-16*%k1
(%i36) r,x=2.0;
(%o36) 2.1510571102112408e-15*%k1
(%i37) r,x=3.0;
(%o37) -(1.1102230246251565e-16*%k2)-5.728750807065808e-14*%k1

[David Billinghurst]

Arfken, Weber and Harris, Mathematical Methods for Physicists 7th ed (2012) Section 14.5 covers the modified Bessel equation pretty well. I see the last line of function modifiedbessel2 should be

return(y = %k1*bessel_i(nu,x-pt) + %k2*bessel_i(-nu,x-pt))))$
This is closer to the bessel2 function and eliminates the issue with ode_check - probably because it is correct.

I will commit this in a day or two.

[David Billinghurst]

Feature added with commit [106905]

[Evert Provoost]

Thank you David! I was able to try it out and that works nicely.

As a side note: I'm alnost done with my experiment and I'll post a small report to the mailing list as suggested earlier.