MATLAB: Problem with symbolic math integration of besselk times exponent

integrationMATLABsymbolicSymbolic Math Toolbox

According to Mathematica and the G&R book, the integral

where K0 is the modified bessel function of the second kind, has a closed form solution that can be expressed using the Meijer G function or the hypergeometric function. However, Matlab's symbolic integrations outputs,

syms x
>> I = int(besselk(0,x/a) * exp(-x/a),x)
 
I =
 
x*exp(-x/a)*(besselk(0, x/a) - besselk(1, x/a))

which is obviously the wrong answer. Matlab itself confirms that, by the way, since,

diff(I,x)
 
ans =
 
exp(-x/a)*(besselk(0, x/a) - besselk(1, x/a)) + x*exp(-x/a)*(besselk(0, x/a)/a - besselk(1, x/a)/a + besselk(1, x/a)/x) - (x*exp(-x/a)*(besselk(0, x/a) - besselk(1, x/a)))/a

What am I doing wrong?

Thanks!

Best Answer

MATLAB is correct, that is one of the ways of expressing the integral.

convert(BesselK(v, z), MeijerG);
         1        /          [[1      1  ]    ]  1  2\
         - MeijerG|[[], []], [[- v, - - v], []], - z |
         2        \          [[2      2  ]    ]  4   /

What am I doing wrong?

You failed to simplify() the differentiation.

Related Solutions

MATLAB: Solving system of equations with parameters

Like sir Walter says I think you should be using ode solvers to solve the problem because it clearly is a system of odes:

tspan=[0 1];  % time span 
ic=[0;1;0;0;0;0;0;0;0;0;0;0]; % initial conditions
[t,x]=ode45(@myod,tspan,ic);  % function call
plot(t,x(:,1)) % likewise another 11 solutions can be plotted
function dSdt = myod(t,x)  % function definition
 S0=x(1);
 S1=x(2);
 S2=x(3);
 S3=x(4);
 S4=x(5);
 S5=x(6);
 S6=x(7);
 S7=x(8);
 S8=x(9);
 S9=x(10);
 S10=x(11);
 S11=x(12);
Ca_sr_conc  = 1.3e-3;
MgATP_conc  = 5e-3;  
MgADP_conc  = 36e-6; 
Pi_conc     = 1e-3;  
Ca=1e-3;  % change it to your value
Ca_cyt_conc = 140e-9; 
k_S0_S1   = 4e7;   
k_S1_S0   = 450;   
k_S1_S2   = 120;   
k_S2_S1   = 25;    
k_S2_S3   = 1e8;   
k_S3_S2   = 16;    
k_S3_S4   = 2e7;  
k_S4_S3   = 40;   
k_S4_S5   = 135;   
k_S5_S4   = 765;   
k_S5_S6   = 15;    
k_S6_S5   = 7.5e4;
k_S6_S7   = 1.5; k_S7_S6   = 2.0;
k_S7_S8   = 350; 
k_S8_S7   = 1.2e5;
k_S8_S9   = 160;   
k_S9_S8   = 250;   
k_S9_S10  = 300;
k_S10_S9  = 9e4;   
k_S10_S11 = 60;    
k_S11_S10 = 60;    
k_S11_S0  = 250;   
k_S0_S11  = 2.7e5; 
dSdt  = [ ((k_S11_S0 * S11) + (k_S1_S0 * S1)) - ((Ca_cyt_conc* k_S0_S1 + Pi_conc * k_S0_S11) * S0);
((Ca_cyt_conc* k_S0_S1 * S0) + (k_S2_S1 * S2)) - ((k_S1_S0 + k_S1_S2) * S1);
((k_S1_S2 * S1) + (k_S3_S2 * S3)) - ((k_S2_S1 + Ca* k_S2_S3) * S2);
((Ca_cyt_conc* k_S2_S3 * S2) + (k_S4_S3 * S4)) - ((k_S3_S2 + MgATP_conc * k_S3_S4) * S3);
((MgATP_conc * k_S3_S4 * S3) + (k_S5_S4 * S5)) - ((k_S4_S3 + k_S4_S5) * S4);
((k_S4_S5 * S4) + (k_S6_S5 * MgADP_conc * S6)) - ((k_S5_S4 + k_S5_S6) * S5);
((k_S5_S6 * S5) + (k_S7_S6 * S7)) - ((k_S6_S5 * MgADP_conc + k_S6_S7) * S6);
((k_S6_S7 * S6) + (Ca_sr_conc * k_S8_S7 * S8)) - ((k_S7_S6 * MgADP_conc + k_S7_S8) * S7);
((k_S7_S8 * S7) + (k_S9_S8 * S9)) - ((Ca_sr_conc * k_S8_S7 + k_S8_S9) * S8);
((k_S8_S9 * S9) + (Ca_sr_conc * k_S10_S9 * S10)) - ((k_S9_S8 + k_S9_S10) * S9);
((k_S9_S10 * S9) + (k_S11_S10 * S11)) - ((k_S10_S11 + Ca_sr_conc * k_S10_S9) * S10);
((k_S10_S11 * S10) + (Pi_conc * k_S0_S11 * S0)) - ((k_S11_S0 + k_S11_S10) * S11)];
end

MATLAB: Negative parameters of Mejier G function

As I indicated before, 0 from the first argument, minus negative 1 from the second argument, gives +1,a positive integer. MeijerG is not defined when one of those differences is a positive integer. There is no way to make such a calculation. You should look up the definition of meijerg and you will see this described.

It is not a limitation of MATLAB or Mathematica or Maple, not a case where "it is just really tough to calculate so we have not bothered": it is fundamental to MeijerG's definition.

Note that it is not a difficulty with negative parameters themselves. You would have the same problem if you for example had 4 in the first parameter and 2 in the second parameter: the difference would be a positive integer and that is not permitted for the function to have meaning. If the first parameter was, for example, 1/2 then -1 would be fine in the second parameter.

Best Answer

Related Solutions

MATLAB: Solving system of equations with parameters

MATLAB: Negative parameters of Mejier G function

Related Question