请教数值积分（已经阅读精华区）

dhp

我想做一个数值积分，数学表达式见我的附件中大括号的积分，被积函数包含两个参数，对其中之一积分，积分上下限是另外一个参数的函数。

我的目标是得到E 相对于 ξ的数值解。

我自己的做法是： 给定ξ的一串值，视η为符号，求出对应的被积函数、上下限，使用int进行积分。运行后有warning，这个查过资料，结果仍然可信。但对于特定的ξ（比如

ξ=1.36）值，会出现错误：The following error occurred converting from sym to double。

大家看看我的思路有没有问题？提提意见吧:因为这个和常见的数值积分不太一样，1是被积函数负责，2是积分上下限是个
ξ的函数。

thx

dhp

我是先考虑得到这个积分值，在考虑整个公式。
clear all;
clc;

syms eta ;

Xi = 1:0.01:10;

Ks = 0.231; Ka = 0.746; Ke = Ks+Ka; g = 0.7; % g is mie scattering parameter K: coefficient
beta_tx = 0; beta_rx = 45*2*pi/360; %transmitter receiver apex angle
theta_tx = pi; theta_rx = 15*2*pi/360; %transmitter receiver half-field view
omega_tx = 4*pi*sin(theta_tx/2); %
theta_s = acos((2-Xi.*Xi-eta*eta)./(Xi.*Xi-eta*eta)); % scattering angle
psi_1 = acos((1+Xi*eta)./(Xi+eta));

phi_r = acos((cos(theta_rx)-cos(beta_rx)*(1+Xi*eta./(Xi+eta)))./(sin(beta_rx)*((Xi.^2-1)*(1-eta^2)).^0.5./(Xi+eta)));
phi_t = acos((cos(theta_tx)-cos(beta_tx)*(1-Xi*eta./(Xi-eta)))./(sin(beta_tx)*((Xi.^2-1)*(1-eta^2)).^0.5./(Xi-eta)));

eta1_r = (Xi*cos(beta_rx+theta_rx)-1)./(Xi-cos(beta_rx+theta_rx));
eta1_t = (1-Xi*cos(beta_tx-theta_tx))./(Xi-cos(beta_tx-theta_tx));
eta2_r = (Xi*cos(beta_rx-theta_rx)-1)./(Xi-cos(beta_rx-theta_rx));
eta2_t = (1-Xi*cos(theta_tx+beta_tx))./(Xi-cos(theta_tx+beta_tx));


g_r = phi_r*cos(beta_rx).*cos(psi_1)+sin(beta_rx)*sin(psi_1).*sin(phi_r); % for isotropic transmitter
g_t = phi_t*cos(beta_rx).*cos(psi_1)+sin(beta_rx)*sin(psi_1).*sin(phi_t); % for isotropic receiver

p_ray = 3*(1+cos(theta_s).*cos(theta_s))/4; % rayleigh scattering phase functin
p_mie = (1-g*g)./(1+g*g-2*g.*cos(theta_s)).^(3/2); % Mie phase function

G_r = 2*g_r.*p_ray./(Xi.*Xi-eta*eta);
G_t = 2*g_t.*p_ray./(Xi.*Xi-eta*eta);
%G_r = 2*g_r.*p_mie./(Xi.*Xi-eta*eta);% for isotropic transmitter
%G_t = 2*g_t.*p_mie./(Xi.*Xi-eta*eta);% for isotropic receiver

E = zeros(1,901);
I = zeros(1,901);

        for n=1:901
E(n) = int(G_r(n), 'eta',eta1_r(n),eta2_r(n));
end;

sigma665

先最外层用一循环,对ξ,给定一范围

这样,积分的上下限也就定了,
然后在对η进行数值积分

dhp

谢谢。我的一个疑问是：我的程序其实也算是一种数值积分了，只不过使用的int命令。您指的数值积分就是用quadl 或quad吧？

还有的是被积函数的表示： 这个被积函数中的每一项都与ξ有关，对于这个函数的定义我感觉有点难度，能指点一二吗？
thx

花如月

int是符号积分吧？

dhp

对，int是符号积分。不过，它也能算出数值，虽然有warning。我对int和quad的区别理解还不深。   
thx

[ 本帖最后由 dhp 于 2008-1-17 09:37 编辑 ]

sigma665

原帖由 dhp 于 2008-1-17 09:29 发表

                               
登录/注册后可看大图

还有的是被积函数的表示： 这个被积函数中的每一项都与ξ有关，对于这个函数的定义我感觉 ...

外循环已经给定了ξ的值,这个已经没有关系了啊

sigma665

原帖由 dhp 于 2008-1-17 09:35 发表

                               
登录/注册后可看大图

对，int是符号积分。不过，它也能算出数值，虽然有warning。我对int和quad的区别理解还不深。   
thx

int算出来的数值是符号型的,如果表达式很复杂,结果也应该是一长串,需要转化成double型

quad 等是数值积分，根据simpson,gauss等数值求积的方法来计算,可以参考数值积分的书

dhp

确实是有关系了，因为被积函数中每项都很长、复杂，所以望而生畏。在我的程序里，其实我是没有定义外部、内部函数的，只是得到一个表达式，循环int积分的。

我试试quad的数值积分吧，有问题继续请教了。

[ 本帖最后由 sigma665 于 2008-1-17 10:06 编辑 ]

dhp

目前我完成了整个数值积分过程，与刚开始的int命令比较，计算结果很快就出来了。不过有warning：Warning: Minimum step size reached; singularity possible.。这个warning后得到的结果，可信不？

我的被积函数是不可积的，非负的。

存在问题：刚开始几个积分结果是复数，这个如何理解？我的被积函数应该是非负的。谢谢。

附我的程序：

function G=myfun_inted(eta,Xi)
%Xi=0.1;
%geometry parameter: angle
g=0.7;

beta_tx = 45*2*pi/360; beta_rx = 45*2*pi/360; %transmitter receiver apex angle
theta_tx = 15*2*pi/360; theta_rx = 15*2*pi/360; %transmitter receiver half-field view

psi_1 = acos((1+Xi.*eta)./(Xi+eta));
theta_s = acos((2-Xi.*Xi-eta.*eta)./(Xi.*Xi-eta.*eta)); % scattering angle



phi_r = acos((cos(theta_rx)-cos(beta_rx)*(1+Xi.*eta./(Xi+eta)))./(sin(beta_rx)*((Xi.^2-1)*(1-eta.^2)).^0.5./(Xi+eta)));
phi_t = acos((cos(theta_tx)-cos(beta_tx)*(1-Xi.*eta./(Xi-eta)))./(sin(beta_tx)*((Xi.^2-1)*(1-eta.^2)).^0.5./(Xi-eta)));

g_r = phi_r*cos(beta_rx).*cos(psi_1)+sin(beta_rx)*sin(psi_1).*sin(phi_r); % for isotropic transmitter
g_t = phi_t*cos(beta_rx).*cos(psi_1)+sin(beta_rx)*sin(psi_1).*sin(phi_t); % for isotropic receiver

p_ray = 3*(1+cos(theta_s).*cos(theta_s))/4; % rayleigh scattering phase functin
p_mie = (1-g*g)./(1+g*g-2*g.*cos(theta_s)).^(3/2); % Mie phase function
G= 2*g_r.*p_ray./(Xi.*Xi-eta.*eta);


% G_r = 2*g_r.*p_ray./(Xi.*Xi-eta*eta);
% G_t = 2*g_t.*p_ray./(Xi.*Xi-eta*eta);

主程序：
clear all;
clc;


beta_tx = 45*2*pi/360; beta_rx = 45*2*pi/360; %transmitter receiver apex angle
theta_tx = 15*2*pi/360; theta_rx = 15*2*pi/360; %transmitter receiver half-field view

%Xi=linspace(1,50);
Xi=1:0.01:10;
for i=1:length(Xi)
   %XX=Xi(i);
eta1_r = (Xi*cos(beta_rx+theta_rx)-1)./(Xi-cos(beta_rx+theta_rx));
eta1_t = (1-Xi*cos(beta_tx-theta_tx))./(Xi-cos(beta_tx-theta_tx));
eta2_r = (Xi*cos(beta_rx-theta_rx)-1)./(Xi-cos(beta_rx-theta_rx));
eta2_t = (1-Xi*cos(theta_tx+beta_tx))./(Xi-cos(theta_tx+beta_tx));
E(i)=quadl(@(eta)myfun_inted(eta,Xi(i)),eta1_r(i),eta2_r(i));
end

[ 本帖最后由 dhp 于 2008-1-18 04:56 编辑 ]

sigma665

singularity possible

积分表达式可能是奇异的...