how can solve a numerical integral?

1 visualización (últimos 30 días)
zara Aghbo
zara Aghbo el 3 de Sept. de 2016
Comentada: Walter Roberson el 3 de Sept. de 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 comentario
Walter Roberson
Walter Roberson el 3 de Sept. de 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.

Iniciar sesión para comentar.

Respuestas (0)

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by