how can solve a numerical integral?
Afficher commentaires plus anciens
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 commentaire
Walter Roberson
le 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.
Réponses (0)
Catégories
En savoir plus sur Symbolic Math Toolbox dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!