how can solve a numerical integral?

7 ビュー (過去 30 日間)
zara Aghbo
zara Aghbo 2016 年 9 月 3 日
コメント済み: Walter Roberson 2016 年 9 月 3 日
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 件のコメント
Walter Roberson
Walter Roberson 2016 年 9 月 3 日
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.

サインインしてコメントする。

回答 (0 件)

カテゴリ

Help Center および File ExchangeCalculus についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by