0
$\begingroup$

I am trying to learn a technique that an author uses for a boundary value problem and, so, I need to be able to use Mathematica to be able to solve the boundary conditions and then plot the final solution so that I can check alongside the journal article. The problem is that the equations involve multiple Bessel functions and the necessity to integrate them. For example, there are four unknown functions similar to the form:

$G(r) = \sum_{n=1}^{\infty} A_{i}h_{j}(\mu_{i}, r) + A_{-i}h_{j}(-\mu_{i}, r)$

where the $h_{j}$ will be linear combinations of the zero order Bessel function and the first order Bessel function and there are a total of 6 different $h_{j}$'s. Then, you have to take a vector of these functions and take the scalar product of it with a vector of the $h_{j}$'s and then finally integrate it over a radius to be able to find the unknown coefficients $A_{i}$ and $A_{-i}$. I attached a couple of screen shots of the relevant pages of one of the journal articles. However, my biggest problem in trying to duplicate the author's results are in using Mathematica to do the final integrals for the coefficients. If I let n=1, the functions will integrate nicely and I can then solve for the coefficients, but for n = 2 or greater, Mathematica stops integrating the functions nicely and just doesn't integrate at all. Then, if I try to NIntegrate, it gives me an error due to the coefficients being unknown. So, I was mostly looking to see if anyone had any advice in potentially making Mathematica able to integrate the equations fully, or, maybe somehow have NIntegrate factor out the coefficients from their grouped terms and then numerically integrate the result, so that I can still get a system of coefficient equations to then solve? I can't find anything related to the latter possibility and I figured that if Mathematica isn't solving the integrals analytically, then it may not be possible to find them analytically, forcing me to need to use NIntegrate. Any help would be appreciated.

enter image description here

enter image description here

  g11[x_, r_] = -BesselJ[1, x * r];

g21[x_, r_] = BesselJ[0, x*r];

g31[x_, r_] =  (c12 - c11) *x* r*BesselJ[1, x*r];

g41[x_, r_] = (c11 - c12) *x* r*BesselJ[0, x*r];

g12[x_, r_] = (3*c11 - c12)/(c11 + c12) *BesselJ[1, x*r] - 
   x * r*BesselJ[0, x*r];

g22[x_, r_] = BesselJ[0, x*r] - x* r*BesselJ[1, x*r];

g32[x_, r_] = (c11 - c12) *x*
   r*((c11 - c12)/(c11 + c12) *BesselJ[1, x*r] - 
     x * r*BesselJ[0, x*r]);

g42[x_, r_] = (c11 - c12) *x*
   r*((c11 + 3*c12)/(c11 + c12) *BesselJ[0, x*r] - 
     x* r*BesselJ[1, x*r]);

U1[x_] = -g31[x, 1]/g32[x, 1];

h1[x_, r_] = g11[x, r] + U1[x]*g12[x, r];

h2[x_, r_] = g21[x, r] + U1[x]*g22[x, r];

h3[x_, r_] = g31[x, r] + U1[x]*g32[x, r];

h4[x_, r_] = g41[x, r] + U1[x]*g42[x, r];

g11[x_, r_] = -BesselJ[1, x * r];

g21[x_, r_] = BesselJ[0, x*r];

g31[x_, r_] =  (c12 - c11) *x* r*BesselJ[1, x*r];

g41[x_, r_] = (c11 - c12) *x* r*BesselJ[0, x*r];

g12[x_, r_] = (3*c11 - c12)/(c11 + c12) *BesselJ[1, x*r] - 
   x * r*BesselJ[0, x*r];

g22[x_, r_] = BesselJ[0, x*r] - x* r*BesselJ[1, x*r];

g32[x_, r_] = (c11 - c12) *x*
   r*((c11 - c12)/(c11 + c12) *BesselJ[1, x*r] - 
     x * r*BesselJ[0, x*r]);

g42[x_, r_] = (c11 - c12) *x*
   r*((c11 + 3*c12)/(c11 + c12) *BesselJ[0, x*r] - 
     x* r*BesselJ[1, x*r]);

U1[x_] = -g31[x, 1]/g32[x, 1];

h1[x_, r_] = g11[x, r] + U1[x]*g12[x, r];

h2[x_, r_] = g21[x, r] + U1[x]*g22[x, r];

h3[x_, r_] = g31[x, r] + U1[x]*g32[x, r];

h4[x_, r_] = g41[x, r] + U1[x]*g42[x, r];

Orz[r_] =  \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m = 
    1\), \(1\)]\((Q[m]*h3[u[\([m]\)], r]\  + \ 
    R[m]*h3[\(-u[\([m]\)]\), r])\)\)

Ozz[r_] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m = 
      1\), \(1\)]\((Q[m]*h4[u[\([m]\)], r]\  + \ 
      R[m]*h4[\(-u[\([m]\)]\), r])\)\) + B0*\[Alpha]*r;

ur[r_] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m = 1\), \(1\)]\((Q[m]*
\*SuperscriptBox[\(E\), \(u[\([m]\)]*l\)]*h1[u[\([m]\)], r]\  + \ R[m]*
\*SuperscriptBox[\(E\), \(\(-u[\([m]\)]\)*l\)]*
       h1[\(-u[\([m]\)]\), r])\)\) + B0*\[Beta]*r;

uz[r_] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m = 1\), \(1\)]\((Q[m]*
\*SuperscriptBox[\(E\), \(u[\([m]\)]*l\)]*h2[u[\([m]\)], r]\  + \ R[m]*
\*SuperscriptBox[\(E\), \(\(-u[\([m]\)]\)*l\)]*
       h2[\(-u[\([m]\)]\), r])\)\) + A0 + B0*l;

F0[r_] = {{0}, {0}, {Orz[r]}, {Ozz[r]}};

F1[r_] = {{ur[r]}, {uz[r]}, {0}, {P}};

Psi1T[m_, r_] = {{h1[x, r], h2[x, r], h3[x, r], h4[x, r]}};

Psi2T[m_, r_] = {{h1[-x, r], h2[-x, r], h3[-x, r], h4[-x, r]}};

Psi0T = {{0, 1, 0, 0}};

Psi01T[r_] = {{\[Beta]*r, 0, 0, \[Alpha]*r}};

(J = {{0, 0, 1, 0}, {0, 0, 0, 1}, {-1, 0, 0, 0}, {0, -1, 0, 
     0}}) // MatrixForm

YY[x_, r_] = Psi2T[x, r].J.(F0[r] - E^(-x*l)*F1[r]);

XX[x_, r_] = Psi1T[x, r].J.(F0[r] - E^(x*l)*F1[r]);

DD[r_] = Psi0T.J.(F0[r] - F1[r]);

UUU[r_] = Psi01T[r].J.(F0[r] - F1[r]);

VVV[r_] = Psi0T.J.F0[r];

EE[r_] = UUU[r] - l*VVV[r];

P = 10;

l = 5;

Coeff1 = Integrate[YY[-u[[1]], r], {r, 0, 1}]
Improve this question
$\endgroup$
9
  • 1
    $\begingroup$ Nothing can be said without the code you tried. So, please help us to help you and include the Mathematica code in copyable form in your post. Moreover, please be more specific about what you need to do and where you are hanging. Quite certainly, nobody will do the whole work for you by reading the article and by programming everything for you. $\endgroup$
    – Henrik Schumacher
    Jul 14, 2018 at 19:32
  • 3
    $\begingroup$ As a side note: Mathematica's Bessel functions (and related functions such as the Hankel functions) are slower by an order of magnitude as they could be. But I doubt that this is the central problem of this post. $\endgroup$
    – Henrik Schumacher
    Jul 14, 2018 at 19:35
  • $\begingroup$ Should I just add the code as a file or something? Also, there wasn't much coding done I'd say, just using Mathematica's build in functions. However, I can post it, if needed. $\endgroup$
    – user571688
    Jul 14, 2018 at 19:36
  • $\begingroup$ Just copy the code to you post, highlight it with the cursor, and press the but "{}" at the top of the editor. $\endgroup$
    – Henrik Schumacher
    Jul 14, 2018 at 19:39
  • $\begingroup$ Okay, I added as much as I thought was relevant. Let me know if you need more. But, the problem is that if I let m = 2 or more for Orz, Ozz, ur, and uz, Integrate[YY] won't fully integrate. $\endgroup$
    – user571688
    Jul 14, 2018 at 19:46

1 Answer 1

Reset to default
1
$\begingroup$

A list u[[m]] is not defined, it is necessary to define it at least as {1,...}. Below is given an example for k = 3.

 u = Table[1, {i, 1, 10}]; k = 3;
g11[x_, r_] = -BesselJ[1, x * r];

g21[x_, r_] = BesselJ[0, x*r];

g31[x_, r_] =  (c12 - c11) *x* r*BesselJ[1, x*r];

g41[x_, r_] = (c11 - c12) *x* r*BesselJ[0, x*r];

g12[x_, r_] = (3*c11 - c12)/(c11 + c12) *BesselJ[1, x*r] - 
      x * r*BesselJ[0, x*r];

g22[x_, r_] = BesselJ[0, x*r] - x* r*BesselJ[1, x*r];

g32[x_, r_] = (c11 - c12) *x*
      r*((c11 - c12)/(c11 + c12) *BesselJ[1, x*r] - 
          x * r*BesselJ[0, x*r]);

g42[x_, r_] = (c11 - c12) *x*
      r*((c11 + 3*c12)/(c11 + c12) *BesselJ[0, x*r] - 
          x* r*BesselJ[1, x*r]);

U1[x_] = -g31[x, 1]/g32[x, 1];

h1[x_, r_] = g11[x, r] + U1[x]*g12[x, r];

h2[x_, r_] = g21[x, r] + U1[x]*g22[x, r];

h3[x_, r_] = g31[x, r] + U1[x]*g32[x, r];

h4[x_, r_] = g41[x, r] + U1[x]*g42[x, r];

g11[x_, r_] = -BesselJ[1, x * r];

g21[x_, r_] = BesselJ[0, x*r];

g31[x_, r_] =  (c12 - c11) *x* r*BesselJ[1, x*r];

g41[x_, r_] = (c11 - c12) *x* r*BesselJ[0, x*r];

g12[x_, r_] = (3*c11 - c12)/(c11 + c12) *BesselJ[1, x*r] - 
      x * r*BesselJ[0, x*r];

g22[x_, r_] = BesselJ[0, x*r] - x* r*BesselJ[1, x*r];

g32[x_, r_] = (c11 - c12) *x*
      r*((c11 - c12)/(c11 + c12) *BesselJ[1, x*r] - 
          x * r*BesselJ[0, x*r]);

g42[x_, r_] = (c11 - c12) *x*
      r*((c11 + 3*c12)/(c11 + c12) *BesselJ[0, x*r] - 
          x* r*BesselJ[1, x*r]);

U1[x_] = -g31[x, 1]/g32[x, 1];

h1[x_, r_] = g11[x, r] + U1[x]*g12[x, r];

h2[x_, r_] = g21[x, r] + U1[x]*g22[x, r];

h3[x_, r_] = g31[x, r] + U1[x]*g32[x, r];

h4[x_, r_] = g41[x, r] + U1[x]*g42[x, r];

Orz[r_] =  \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m = 
     1\), \(k\)]\((Q[m]*h3[u[\([\)\(m\)\(]\)], r]\  + \ 
     R[m]*h3[\(-u[\([\)\(m\)\(]\)]\), r])\)\);

Ozz[r_] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m = 
      1\), \(k\)]\((Q[m]*h4[u[\([\)\(m\)\(]\)], r]\  + \ 
      R[m]*h4[\(-u[\([\)\(m\)\(]\)]\), r])\)\) + B0*\[Alpha]*r;

ur[r_] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m = 1\), \(k\)]\((Q[m]*
\*SuperscriptBox[\(E\), \(u[\([\)\(m\)\(]\)]*l\)]*
       h1[u[\([\)\(m\)\(]\)], r]\  + \ R[m]*
\*SuperscriptBox[\(E\), \(\(-u[\([\)\(m\)\(]\)]\)*l\)]*
       h1[\(-u[\([\)\(m\)\(]\)]\), r])\)\) + B0*\[Beta]*r;

uz[r_] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(m = 1\), \(k\)]\((Q[m]*
\*SuperscriptBox[\(E\), \(u[\([\)\(m\)\(]\)]*l\)]*
       h2[u[\([\)\(m\)\(]\)], r]\  + \ R[m]*
\*SuperscriptBox[\(E\), \(\(-u[\([\)\(m\)\(]\)]\)*l\)]*
       h2[\(-u[\([\)\(m\)\(]\)]\), r])\)\) + A0 + B0*l;

F0[r_] = {{0}, {0}, {Orz[r]}, {Ozz[r]}};

F1[r_] = {{ur[r]}, {uz[r]}, {0}, {P}};

Psi1T[m_, r_] = {{h1[x, r], h2[x, r], h3[x, r], h4[x, r]}};

Psi2T[m_, r_] = {{h1[-x, r], h2[-x, r], h3[-x, r], h4[-x, r]}};

Psi0T = {{0, 1, 0, 0}};

Psi01T[r_] = {{\[Beta]*r, 0, 0, \[Alpha]*r}};

(J = {{0, 0, 1, 0}, {0, 0, 0, 1}, {-1, 0, 0, 0}, {0, -1, 0, 
          0}});

YY[x_, r_] = Psi2T[x, r].J.(F0[r] - E^(-x*l)*F1[r]);

XX[x_, r_] = Psi1T[x, r].J.(F0[r] - E^(x*l)*F1[r]);

DD[r_] = Psi0T.J.(F0[r] - F1[r]);

UUU[r_] = Psi01T[r].J.(F0[r] - F1[r]);

VVV[r_] = Psi0T.J.F0[r];

EE[r_] = UUU[r] - l*VVV[r];

P = 10;

l = 5;

Coeff1 = Integrate[YY[-u[[1]], r], {r, 0, 1}] // AbsoluteTiming
Out[47]: {95.9164, {{(-2 (c11^2 - c12^2) (-3 (c11 + c12) BesselJ[0, 1]^4 + 
         6 c12 BesselJ[0, 1]^3 BesselJ[1, 1] + 
         4 (2 c11 - c12) BesselJ[0, 1]^2 BesselJ[1, 1]^2 - 
         2 (c11 - 2 c12) BesselJ[0, 1] BesselJ[1, 
           1]^3 + (7 c11 + c12) BesselJ[1, 1]^4) (Q[1] + Q[2] + Q[3] +
          R[1] + R[2] + R[3]) + 
      2 (c11 - c12) (c11 + c12) (3 (c11 + c12) BesselJ[0, 1]^4 - 
         6 (2 c11 + c12) BesselJ[0, 1]^3 BesselJ[1, 1] + 
         4 (4 c11 + c12) BesselJ[0, 1]^2 BesselJ[1, 1]^2 - 
         2 (5 c11 + 2 c12) BesselJ[0, 1] BesselJ[1, 
           1]^3 + (5 c11 - c12) BesselJ[1, 
           1]^4) (E^10 (Q[1] + Q[2] + Q[3]) + R[1] + R[2] + R[3]) + 
      3 (c11^2 - 
         c12^2) (12 B0 (c11 - c12) E^5 \[Beta] BesselJ[1, 1]^3 - 
         BesselJ[0, 
           1]^2 (-B0 (c11 + c12) E^5 \[Beta] (8 BesselJ[1, 1] - 
               Hypergeometric0F1Regularized[3, -(1/4)]) + 
            2 (3 c11 - c12) BesselJ[1, 
              1]^2 (E^10 (Q[1] + Q[2] + Q[3]) - R[1] - R[2] - R[3])) +
          BesselJ[0, 1] BesselJ[1, 
           1] (-B0 E^5 \[Beta] (20 c11 BesselJ[1, 1] + 
               4 c12 BesselJ[1, 1] - 
               c11 Hypergeometric0F1Regularized[3, -(1/4)] + 
               c12 Hypergeometric0F1Regularized[
                 3, -(1/4)]) + (4 c11 BesselJ[1, 1]^2 - 
               c11 HypergeometricPFQ[{3/2}, {2, 3}, -1] - 
               c12 HypergeometricPFQ[{3/2}, {2, 3}, -1]) (E^10 (Q[1] +
                   Q[2] + Q[3]) - R[1] - R[2] - R[3])) + 
         2 (c11 + c12) BesselJ[0, 
           1]^4 (E^10 (Q[1] + Q[2] + Q[3]) - R[1] - R[2] - R[3]) - 
         BesselJ[1, 
           1]^2 (c11 (HypergeometricPFQ[{3/2}, {2, 3}, -1] - 
               HypergeometricPFQ[{1/2, 2}, {1, 1, 3}, -1]) - 
            c12 HypergeometricPFQ[{1/2, 2}, {1, 1, 
               3}, -1]) (E^10 (Q[1] + Q[2] + Q[3]) - R[1] - R[2] - 
            R[3]) + 
         4 c12 BesselJ[0, 1]^3 BesselJ[1, 
           1] (-E^10 (Q[1] + Q[2] + Q[3]) + R[1] + R[2] + R[3])) + 
      2 (3 (c11 - c12) (c11 + c12)^2 BesselJ[0, 
           1]^4 (Q[1] + Q[2] + Q[3] - R[1] - R[2] - R[3]) + (5 c11^3 -
             c11^2 c12 - 5 c11 c12^2 + c12^3) BesselJ[1, 
           1]^4 (Q[1] + Q[2] + Q[3] - R[1] - R[2] - R[3]) - 
         2 (c11 - c12) BesselJ[1, 
           1]^3 (6 B0 c12 \[Alpha] + 
            5 c11^2 BesselJ[0, 
              1] (Q[1] + Q[2] + Q[3] - R[1] - R[2] - R[3]) + 
            c12 (2 c12 BesselJ[0, 
                 1] (Q[1] + Q[2] + Q[3] - R[1] - R[2] - R[3]) - 
               15 E^5 \[Pi] StruveH[0, 1]) + 
            c11 (7 c12 BesselJ[0, 
                 1] (Q[1] + Q[2] + Q[3] - R[1] - R[2] - R[3]) + 
               15 E^5 \[Pi] StruveH[0, 1])) + 
         30 (c11 + c12)^2 E^5 BesselJ[0, 
           1]^3 (-2 + \[Pi] StruveH[1, 1]) - 
         6 (c11 + c12) BesselJ[0, 1]^2 BesselJ[1, 
           1] (2 c11^2 BesselJ[0, 
              1] (Q[1] + Q[2] + Q[3] - R[1] - R[2] - R[3]) + 
            c12 (c12 BesselJ[0, 
                 1] (-Q[1] - Q[2] - Q[3] + R[1] + R[2] + R[3]) + 
               5 E^5 (2 + \[Pi] StruveH[0, 1] - 
                  2 \[Pi] StruveH[1, 1])) + 
            c11 (c12 BesselJ[0, 
                 1] (-Q[1] - Q[2] - Q[3] + R[1] + R[2] + R[3]) + 
               5 E^5 (-6 + \[Pi] StruveH[0, 1] + 
                  2 \[Pi] StruveH[1, 1]))) + 
         2 BesselJ[0, 1] BesselJ[1, 
           1]^2 (8 c11^3 BesselJ[0, 
              1] (Q[1] + Q[2] + Q[3] - R[1] - R[2] - R[3]) + 
            c12^2 (6 B0 \[Alpha] - 
               2 c12 BesselJ[0, 
                 1] (Q[1] + Q[2] + Q[3] - R[1] - R[2] - R[3]) + 
               15 E^5 \[Pi] (-2 StruveH[0, 1] + StruveH[1, 1])) + 

            2 c11 c12 (3 B0 \[Alpha] + 
               4 c12 BesselJ[0, 
                 1] (-Q[1] - Q[2] - Q[3] + R[1] + R[2] + R[3]) - 
               15 E^5 (-2 + \[Pi] StruveH[1, 1])) + 
            c11^2 (2 c12 BesselJ[0, 
                 1] (Q[1] + Q[2] + Q[3] - R[1] - R[2] - R[3]) + 
               15 E^5 (-4 + \[Pi] (2 StruveH[0, 1] + 
                    StruveH[1, 1]))))))/(12 ((c11 + c12) BesselJ[0, 
          1] - (c11 - c12) BesselJ[1, 1])^2)}}}
Improve this answer
$\endgroup$
3
  • $\begingroup$ Hi, thank you for the reply. I was able to get it to integrate by your method, however, I noticed that once I set u to something like u = Table[i^2, {i, 1, 10}], so that the elements vary, the equation stopped integrating again. In my actual problem they would be varying. $\endgroup$
    – user571688
    Jul 16, 2018 at 23:39
  • $\begingroup$ Can you publish a link to the article, where are the equations taken from? I can not see the drawings that you published. $\endgroup$
    – Alex Trounev
    Jul 17, 2018 at 4:02
  • $\begingroup$ kundoc.com/… So, I attached a url that I found one of the articles at (I originally got them through my school's library), but the problem in the article isn't quite the same as the one I am trying to figure out. The article I am trying to figure out has no gravity included, but the equivalent equations that I am having difficulty integrating would be 85-88, but for my case the right hand side would be all 0. The relevant equations in the code are YY, XX, DD, and EE. $\endgroup$
    – user571688
    Jul 17, 2018 at 11:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.