1. General
  2. Spørgsmål
0
Besvaret

Definite integrals with mp-objects

Jeremy 9 år siden • opdateret af Pavel Holoborodko 9 år siden • 10

Hello,


Is there a possibility to calculate definite integrals of sum of products of sines/cosines with cosh/sinh with mp-objects?


Kind regards,


Jeremy

Under vurdering
Pavel Holoborodko 9 år siden

Hello Jeremy,


Yes, of course. Toolbox provides full spectrum of numerical integration routines - from fixed Gauss quadrature to adaptive quadgk.


Or did you mean something else?

Jeremy 9 år siden

For example, I would like to calculate the definite integral "int(sin(b*x)*cosh(c*x)" in the variable x over the interval [0,L], in which b,c and L are constant values (which I defined as mp-objects).

Pavel Holoborodko 9 år siden

Define your integrand as a function, call quadgk with corresponding parameters.


Make sure the values 0, L and tolerance are also of mp-type. Toolbox will detect the mp-objects and will call its own multiprecision version of quadgk.


Let me know how it works.

Pavel Holoborodko 9 år siden

Here is simple example for demonstration: 

mp.Digits(50);
f = @(x)sin(x);
[q,errbnd] = quadgk(f,mp(0),mp(1),'RelTol',100*mp.eps,'AbsTol',100*mp.eps, 'MaxIntervalCount', 1000000)
q = 
    0.45969769413186028259906339255702339626768957938206
errbnd = 
    7.872028207137840807482477381844110449433981844857e-51
Just use your own integrand and interval boundaries
Jeremy 9 år siden

Can I also use this to create a matrix of solutions, since I need to calculate these integrals with varying coefficients b and c?

Besvaret
Pavel Holoborodko 9 år siden

Yes, of course, toolbox follows standard Matlab syntax and rules.


If you can do something with standard Matlab - most likely it can be done with toolbox too.

Pavel Holoborodko 9 år siden

Hi again,


I couldn't find your last comment (did you delete it?).


Just few comments on the code.

The function should be able to accept the vector-parameters (quadqk is a vectorized method). So instead of 

f = @(x)(cos(p*x)+cosh(p*x)+q*sin(p*x)+q*sinh(p*x))*(cos(r*x)+cosh(r*x)+s*sin(r*x)+s*sinh(r*x));

use element-wise multiplication: 

f = @(x)(cos(p*x)+cosh(p*x)+q*sin(p*x)+q*sinh(p*x)).*(cos(r*x)+cosh(r*x)+s*sin(r*x)+s*sinh(r*x));

Also be careful with precision and interval counts. My example shows extreme values, please use more reasonable at first, e.g.: 

mp.Digits(34);
...
[q,errbnd] = quadgk(f,mp(0),L,'RelTol',100*mp.eps(),'AbsTol',100*mp.eps(), 'MaxIntervalCount', 650)

Quadgk will warn you if it needs more levels of recursion (MaxIntervalCount)

Jeremy 9 år siden

Thank you for your assistance! You've been of great help!

Jeremy 9 år siden

Hi again,


I can't seem to numerically stabilize my code (see below). I expect the diagonal elements in the matrix "Ortho" to be equal to L (here = 10) and the off-diagonal elements equal to zero. I'm trying to prove the orthogonality of the function f...


mp.Digits(20);

L = mp(10);
N = mp(5);
fun = @(x) cos(x)*cosh(x)-1;
B = mp(zeros(N,1));
B(1,1) = fzero(fun,[4,5]);
B(2,1) = fzero(fun,[7,8]);
for k=3:N
B(k,1) = (2*(k+1)-1)*pi/2;
end
for i=1:N
b(i,1) = B(i,1)/L;
end
ksi = mp(zeros(N,1));
for i=1:N
r = B(i,1);
ksi(i,1) = -(cos(r)-cosh(r))/(sin(r)-sinh(r));
end
Ortho = mp(zeros(N,N));
for i=1:N
p = b(i,1);
q = ksi(i,1);
for j=1:N
r = b(j,1);
s = ksi(j,1);
f = @(x)(cos(p*x)+cosh(p*x)+q*sin(p*x)+q*sinh(p*x)).*(cos(r*x)+cosh(r*x)+s*sin(r*x)+s*sinh(r*x));
[q,errbnd] = quadgk(f,mp(0),L,'RelTol',100*mp.eps(),'AbsTol',100*mp.eps(), 'MaxIntervalCount', 100000);
Ortho(j,i) = q;
end
end



Pavel Holoborodko 9 år siden

I am not sure, probably you need to use better precision (from numerical point of view).

But why don't you compute the integral analytically?


Try using Maple or Mathematica. I have just tried Maple - it gave long but sensible expression for the any r,p,q,s.


Kundesupport af UserEcho