Integral of another integral

I want to integrate h1b_eta and then plot it (m1b_eta). However, h1b_eta depends on other functions that also contain integrals. How can I do this?

h1b_eta=@(eta) blablabla and some integrals;

m1b_eta=@(eta)  pi*(airy(2,eta)*integral(@(n) airy(n)*h1b_eta(n),eta_infts,eta)-airy(eta)*integral(@(n) airy(2,n)*h1b_eta(n),eta0ts,eta));

V1VEC=arrayfun(m1b_eta,0:0.01:N);
plot(abs(V1VEC),0:0.01:N)

I want to do this symbollically because I need to integrate m1b_eta aftwerwards (this can now be done numerically)

[carlos g]

@Walter Roberson How can I use "result"? I am not that handy with symbolic functions. A is basically a variable (a number) and B is basically my variable eta, the same for GAMMA, that you defined

Answer 1

[Walter Roberson]

Code h1b_eta as a symbolic function; then

syms m1b_eta(eta) n
Pi = sym('pi');
m1b_eta(eta) = Pi * (airy(2,eta)*int(airy(n)*h1b_eta(n),n,eta_infts,eta) - airy(eta)*int( airy(2,n)*h1b_eta(n),n,eta0ts,eta));

x = 0:0.01:N;
V1VEC = m1b_eta(x);
plot(x, V1VEC)

Possibly the call will not be enough because it might have to turn the integral into numeric; you might need to use

V1VEC = arrayfun(@(ETA) double(m1b_eta(ETA)), x);

Also, if your expression is complicated, you might wan to use vpaintegral() instead of int() for performance: if you do, be sure to specify the RelTol option to get enough accuracy for your purpose.

[carlos g]

I still get problems, this time with integral2, which is used by some of the functions. I include the whole code here

syms x
syms y
syms h1b_etas(eta) n
syms m1b_etas(eta) n
Pi = sym('pi');

ai2s(x)=x*airy(x);
bi2s(x)=x*airy(2,x);

Gi1s(x)= -(airy(2,x)*int(airy(n),n,eta_inf1,x)-airy(x)*int(airy(2,n),n,eta01,x));
Gip1s(x)= -(-airy(3,x)*int(airy(n),n,x,eta_inf1)-airy(1,x)...
        *int(airy(2,n),n,eta01,x));
Gipp1s(x)= -(bi2s(x)*int(airy(n),n,eta_inf1,x)+airy(3,x)...
        *airy(x)-ai2s(x)*int(airy(2,n),n,eta01,x)...
        -airy(1,x)*airy(2,x));
Gi2s(x)= -(airy(2,x)*int(airy(n),n,eta_inf2,x)-airy(x)*int(airy(2,n),n,eta02,x));
Gip2s(x)= -(-airy(3,x)*int(airy(n),n,x,eta_inf2)-airy(1,x)...
        *int(airy(2,n),n,eta02,x));
Gipp2s(x)= -(bi2s(x)*int(airy(n),n,eta_inf2,x)+airy(3,x)...
        *airy(x)-ai2s(x)*int(airy(2,n),n,eta02,x)...
        -airy(1,x)*airy(2,x));

V1s(eta)=((i*alpha1*lambda_0)^(-1/3))*q1*integral2(@(x, y) airy(y),...
        eta01, eta, eta01, @(x) x);
V2s(eta)=((i*alpha2*lambda_0)^(-1/3))*q2*integral2(@(x, y) airy(y),...
        eta02, eta, eta02, @(x) x);

eta1s(y)= ((i*alpha1*lambda_0)^(1/3))*y+eta01;
eta2s(y)= ((i*alpha2*lambda_0)^(1/3))*y+eta02;
etatss(y)= ((i*alphats*lambda_0)^(1/3))*y+eta0ts;

DV1s(eta)= q1*int(airy(n),n,eta01,eta);
DV2s(eta)= q2*int(airy(n),n,eta02,eta);
W1s(eta)= -i*beta1*airy_D1/(alpha1^2+beta1^2)*Pi*Gi1s(eta);
U1s(eta)= i/alpha1*DV1s(eta)-beta1/alpha1*W1s(eta);
W2s(eta)= -i*beta2*airy_D2/(alpha2^2+beta2^2)*Pi*Gi2s(eta);
U2s(eta)= i/alpha2*DV2s(eta)-beta2/alpha2*W2s(eta);

DDV1s(eta)= ((i*alpha1*lambda_0)^(1/3))*q1*airy(eta);
DDV2s(eta)= ((i*alpha2*lambda_0)^(1/3))*q2*airy(eta);

DDDV1(eta)= ((i*alpha1*lambda_0)^(2/3))*q1*airy(1,eta);
DDDV2(eta)= ((i*alpha2*lambda_0)^(2/3))*q2*airy(1,eta);

DW1s(eta)= -i*beta1*airy_D1*Pi*((i*alpha1*lambda_0)^(1/3))*Gip1s(eta);
DU1s(eta)= i/alpha1*DDV1s(eta)-beta1/alpha1*DW1s(eta);
DDW1s(eta)= -i*beta1*airy_D1*Pi*((i*alpha1*lambda_0)^(2/3))*Gipp1s(eta);
DDU1s(eta)= i/alpha1*DDDV1s(eta)-beta1/alpha1*DDW1s(eta);
DW2s(eta)= -i*beta2*airy_D2*Pi*((i*alpha2*lambda_0)^(1/3))*Gip2s(eta);
DU2s(eta)= i/alpha2*DDV2s(eta)-beta2/alpha2*DW2s(eta);
DDW2s(eta)= -i*beta2*airy_D2*Pi*((i*alpha2*lambda_0)^(2/3))*Gipp2s(eta);
DDU2s(eta)= i/alpha2*DDDV2s(eta)-beta2/alpha2*DDW2s(eta);

N1Y_etas(eta)=i*alphats*(DU1s(eta)*conj(U2s(eta))+U1s(eta)...
    *conj(DU2s(eta)))+i*beta1*(conj(DW2s(eta))*U1s(eta)...
    +conj(W2s(eta))*DU1s(eta))-i*beta2*(DW1s(eta)...
    *conj(U2s(eta))+W1s(eta)*conj(DU2s(eta)))...
    +(DV1s(eta)*conj(DU2s(eta))+V1s(eta)...
    *conj(DDU2s(eta)))+conj(DV2s(eta))*DU1s(eta)...
    +conj(V2s(eta))*DDU1s(eta);

N3Y_etas(eta)=i*alpha1*(DW1s(eta)*conj(U2s(eta))+W1s(eta)...
    *conj(DU2s(eta)))-i*alpha2*(conj(DW2s(eta))*U1s(eta)...
    +conj(W2s(eta))*DU1s(eta))+i*betats*(DW1s(eta)...
    *conj(DW2s(eta))+W1s(eta)*conj(DW2s(eta)))...
    +(DV1s(eta)*conj(DW2s(eta))+V1s(eta)...
    *conj(DDW2s(eta)))+(conj(DV2s(eta))*DW1s(eta)...
    +conj(V2s(eta))*DDW1s(eta));

h1b_etas(eta)= -(1/(lambda_0*alphats))*(alphats*N1Y_etas(eta)+betats*N3Y_etas(eta));

m1b_etas(eta) = Pi * (airy(2,eta)*int(airy(n)*h1b_etas(n),n,eta_infts,eta) - airy(eta)*int( airy(2,n)*h1b_etas(n),n,eta0ts,eta));
vec= 0:0.01:N;
V1VEC = m1b_etas(x);
plot(vec, V1VEC)

[Walter Roberson]

integral2(@(x, y) airy(y), eta01, eta, eta01, @(x) x)

is

int(int(airy(y), y, eta01, x), x, eta01, eta)

You use the same form but different variables in the next line, so it might help to know that the general form

int(int(airy(y),y,A,x),x,A,B)

has a closed-form solution,

   syms A B; GAMMA=@(x) gamma(sym(x)); Pi = sym('pi');
   result = (1/sym(120))*(20*A*Pi*sym(3)^sym(1/3)*(A-2*B)*hypergeom(sym([1/3]), sym([4/3, 5/3]), (1/sym(9))*A^3)-30*A^2*sym(3)^(1/sym(6))*GAMMA(2/3)^2*(A-B)*hypergeom(sym([2/3]), sym([4/3, 5/3]), (1/sym(9))*A^3)+A^4*Pi*sym(3)^sym(1/3)*(A-B)*hypergeom(sym([4/3]), sym([7/3, 8/3]), (1/sym(9))*A^3)+20*B^2*sym(3)^sym(1/3)*Pi*hypergeom(sym([1/3]), sym([4/3, 5/3]), (1/9)*B^3)-60*sym(3)^sym(1/6)*GAMMA(2/3)^2*hypergeom(sym([-1/3]), sym([1/3, 2/3]), (1/9)*A^3)+60*sym(3)^sym(1/6)*GAMMA(2/3)^2*hypergeom(sym([-1/3]), sym([1/3, 2/3]), (1/9)*B^3))/(Pi*GAMMA(2/3))

[carlos g]

@Walter Roberson How can I use "result"? I am not that handy with symbolic functions. A is basically a variable (a number) and B is basically my variable eta, the same for GAMMA, that you defined

[Walter Roberson]

V1s(eta)=((i*alpha1*lambda_0)^(-1/3))*q1*integral2(@(x, y) airy(y),...
      eta01, eta, eta01, @(x) x);
V2s(eta)=((i*alpha2*lambda_0)^(-1/3))*q2*integral2(@(x, y) airy(y),...
      eta02, eta, eta02, @(x) x);

comes out as

GAMMA=@(x) gamma(sym(x)); Pi = sym('pi');

V1s(eta) = (1/3)*(eta-eta01)*(-(3/4)*3^(1/6)*(hypergeom([2/3], [4/3, 5/3], (1/9)*x^3)*x^2-eta01^2*hypergeom([2/3], [4/3, 5/3], (1/9)*eta01^3))*GAMMA(2/3)^2+Pi*3^(1/3)*(x*hypergeom([1/3], [2/3, 4/3], (1/9)*x^3)-eta01*hypergeom([1/3], [2/3, 4/3], (1/9)*eta01^3)))*q1/((I*alpha1*lambda_0)^(1/3)*Pi*GAMMA(2/3));
V2s(eta) = (1/6)*(-2*3^(1/3)*eta02*Pi*(eta-(1/2)*eta02)*hypergeom([1/3], [4/3, 5/3], (1/9)*eta02^3)+(3/2)*GAMMA(2/3)^2*3^(1/6)*eta02^2*(eta-eta02)*hypergeom([2/3], [4/3, 5/3], (1/9)*eta02^3)-(1/20)*Pi*3^(1/3)*eta02^4*(eta-eta02)*hypergeom([4/3], [7/3, 8/3], (1/9)*eta02^3)+Pi*hypergeom([1/3], [4/3, 5/3], (1/9)*eta^3)*3^(1/3)*eta^2+3*GAMMA(2/3)^2*hypergeom([-1/3], [1/3, 2/3], (1/9)*eta^3)*3^(1/6)-3*GAMMA(2/3)^2*3^(1/6)*hypergeom([-1/3], [1/3, 2/3], (1/9)*eta02^3))*q2/((I*alpha2*lambda_0)^(1/3)*Pi*GAMMA(2/3));