请教数值积分（已经阅读精华区）
我想做一个数值积分，数学表达式见我的附件中大括号的积分，被积函数包含两个参数，对其中之一积分，积分上下限是另外一个参数的函数。我的目标是得到E 相对于 ξ的数值解。
我自己的做法是： 给定ξ的一串值，视η为符号，求出对应的被积函数、上下限，使用int进行积分。运行后有warning，这个查过资料，结果仍然可信。但对于特定的ξ（比如
ξ=1.36）值，会出现错误：The following error occurred converting from sym to double。
大家看看我的思路有没有问题？提提意见吧:因为这个和常见的数值积分不太一样，1是被积函数负责，2是积分上下限是个
ξ的函数。
thx
附我的matlab程序
我是先考虑得到这个积分值，在考虑整个公式。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;
回复 楼主 的帖子
先最外层用一循环,对ξ,给定一范围这样,积分的上下限也就定了,
然后在对η进行数值积分
回复 3楼 的帖子
谢谢。我的一个疑问是：我的程序其实也算是一种数值积分了，只不过使用的int命令。您指的数值积分就是用quadl 或quad吧？还有的是被积函数的表示： 这个被积函数中的每一项都与ξ有关，对于这个函数的定义我感觉有点难度，能指点一二吗？
thx
回复 4楼 的帖子
int是符号积分吧？回复 5楼 的帖子
对，int是符号积分。不过，它也能算出数值，虽然有warning。我对int和quad的区别理解还不深。thx
[ 本帖最后由 dhp 于 2008-1-17 09:37 编辑 ] 原帖由 dhp 于 2008-1-17 09:29 发表 http://www.chinavib.com/forum/images/common/back.gif
还有的是被积函数的表示： 这个被积函数中的每一项都与ξ有关，对于这个函数的定义我感觉 ...
外循环已经给定了ξ的值,这个已经没有关系了啊 原帖由 dhp 于 2008-1-17 09:35 发表 http://www.chinavib.com/forum/images/common/back.gif
对，int是符号积分。不过，它也能算出数值，虽然有warning。我对int和quad的区别理解还不深。
thx
int算出来的数值是符号型的,如果表达式很复杂,结果也应该是一长串,需要转化成double型
quad 等是数值积分，根据simpson,gauss等数值求积的方法来计算,可以参考数值积分的书
回复 7楼 的帖子
确实是有关系了，因为被积函数中每项都很长、复杂，所以望而生畏。在我的程序里，其实我是没有定义外部、内部函数的，只是得到一个表达式，循环int积分的。我试试quad的数值积分吧，有问题继续请教了。
[ 本帖最后由 sigma665 于 2008-1-17 10:06 编辑 ]
继续昨天的题目请教数值积分！
目前我完成了整个数值积分过程，与刚开始的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 编辑 ]
回复 楼主 的帖子
singularity possible积分表达式可能是奇异的...
页:
[1]