how can solve a numerical integral?

8 views (last 30 days)
zara Aghbo
zara Aghbo on 3 Sep 2016
Commented: Walter Roberson on 3 Sep 2016
I try to solve a numerical fourth order Integral. But it has taken a lot of time without any solution or even warning. Could you please help me to find it's problem. My code is
clc
clear all
close all
syms r R b r1 r2 r3 r4 r5 r6 z z1 z2 z3 z4 z5 z6
%syms r R b r1 r2 z z1 z2
hc=197.326;
alphas=1.54;
a=25.13;
l=-8/3;
s=1;
e=0.5;
M=313;
b=0.603;
f1=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*(r.^2+(s./2).^2-r.*s.*cos(z))));
f2=(1./pi.*b.^2)^(3./4).*(exp((-1./b.^2).*(r.^2+(s./2).^2+r.*s.*cos(z))));
k=(int(int(f1.*f2,'r',0,inf ),'z',0,pi));
N=(1+e.^2+2.*e.*k).^(1./2);
fl_1(r1)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r1).^2+(s./2).^2-(r1).*s.*cos(z1))));
fl_2(r2)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r2).^2+(s./2).^2-(r2).*s.*cos(z2))));
fl_3(r3)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r3).^2+(s./2).^2-(r3).*s.*cos(z3))));
fl_4(r4)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r4).^2+(s./2).^2-(r4).*s.*cos(z4))));
fl_5(r5)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r5).^2+(s./2).^2-(r5).*s.*cos(z5))));
fl_6(r6)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r6).^2+(s./2).^2-(r6).*s.*cos(z6))));
fR_1(r1)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r1).^2+(s./2).^2+(r1).*s.*cos(z1))));
fR_2(r2)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r2).^2+(s./2).^2+(r2).*s.*cos(z2))));
fR_3(r3)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r3).^2+(s./2).^2+(r3).*s.*cos(z3))));
fR_4(r4)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r4).^2+(s./2).^2+(r4).*s.*cos(z4))));
fR_5(r5)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r5).^2+(s./2).^2+(r5).*s.*cos(z5))));
fR_6(r6)=(1./pi.*b.^2).^(3./4).*(exp((-1./b.^2).*((r6).^2+(s./2).^2+(r6).*s.*cos(z6))));
uL_1(r1)=(fl_1(r1)+e.*fR_1(r1))./N;
uL_2(r2)=(fl_2(r2)+e.*fR_2(r2))./N;
uL_3(r3)=(fl_3(r3)+e.*fR_3(r3))./N;
uR_4(r4)=(fR_4(r4)+e.*fl_4(r4))./N;
uR_5(r5)=(fR_5(r5)+e.*fl_5(r5))./N;
uR_6(r6)=(fR_6(r6)+e.*fl_6(r6))./N;
g1=uL_1(r1).*uL_2(r2).*(1./(r1-r2)).*uL_1(r1).*uL_2(r2).*(r1).^2.*(r2).^2.*sin(z1).*sin(z2);
g2= matlabFunction(g1);
G_11= integral2(@(r1,r2)arrayfun(@(r1,r2)integral2(@(z1,z2)(g2(r1,r2,z1,z2)),0,pi,0,pi),r1,r2),1,2,@(r1)r1,2);
G_12= integral2(@(r1,r2)arrayfun(@(r1,r2)integral2(@(z1,z2)(g2(r1,r2,z1,z2)),0,pi,0,pi),r1,r2),1,2,1,@(r1)r1);
G_1=G_11+G_12;
G_2=3.*(l./4).*(hc).*(alphas).*(2.*pi).^2.*(G_1);
g4=uR_4(r4).*uR_6(r6).*(1./(r4-r6)).*uR_4(r4).*uR_6(r6).*(r4).^2.*(r6).^2.*sin(z4).*sin(z6);
g6= matlabFunction(g4);
G_44= integral2(@(r4,r6)arrayfun(@(r4,r6)integral2(@(z4,z6)(g6(r4,r6,z4,z6)),0,pi,0,pi),r4,r6),1,2,@(r4)r4,2);
G_46= integral2(@(r4,r6)arrayfun(@(r4,r6)integral2(@(z4,z6)(g6(r4,r6,z4,z6)),0,pi,0,pi),r4,r6),1,2,1,@(r4)r4);
G_4=G_44+G_46;
G_6=3.*(l./4).*(hc).*(alphas).*(2.*pi).^2.*(G_4);
g3=uL_3(r3).*uR_5(r5).*(1./(r3-r5)).*uL_3(r3).*uR_5(r5).*(r3).^2.*(r5).^2.*sin(z3).*sin(z5);
g5= matlabFunction(g3);
G_33= integral2(@(r3,r5)arrayfun(@(r3,r5)integral2(@(z3,z5)(g5(r3,r5,z3,z5)),0,pi,0,pi),r3,r5),1,2,@(r3)r3,2);
G_35= integral2(@(r3,r5)arrayfun(@(r3,r5)integral2(@(z3,z5)(g5(r3,r5,z3,z5)),0,pi,0,pi),r3,r5),1,2,1,@(r3)r3);
G_3=G_33+G_35;
G_5=3.*(l./4).*(hc).*(alphas).*(2.*pi).^2.*(G_3);
V_g=G_2+G_6+G_5
  1 Comment
Walter Roberson
Walter Roberson on 3 Sep 2016
Is it correct that G_11 is
integral of g1, z1 = 0 .. Pi, z2 = 0 .. Pi, r2 = r1 .. 2, r1 = 1 .. 2)
?
Maple is telling me that the numerator of that is undefined.

Sign in to comment.

Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!