how can solve a numerical integral?
6 vues (au cours des 30 derniers jours)
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)
Voir également
Catégories
En savoir plus sur Calculus dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!