声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1861|回复: 10

[编程技巧] 请教数值积分(已经阅读精华区)

[复制链接]
发表于 2008-1-17 09:16 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

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

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

ξ=1.36)值,会出现错误:The following error occurred converting from sym to double

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

ξ的函数。

thx
算法reilly.jpg
回复
分享到:

使用道具 举报

 楼主| 发表于 2008-1-17 09:20 | 显示全部楼层

附我的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;
发表于 2008-1-17 09:23 | 显示全部楼层

回复 楼主 的帖子

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

这样,积分的上下限也就定了,
然后在对η进行数值积分
 楼主| 发表于 2008-1-17 09:29 | 显示全部楼层

回复 3楼 的帖子

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

还有的是被积函数的表示: 这个被积函数中的每一项都与ξ有关,对于这个函数的定义我感觉有点难度,能指点一二吗?
thx
发表于 2008-1-17 09:34 | 显示全部楼层

回复 4楼 的帖子

int是符号积分吧?
 楼主| 发表于 2008-1-17 09:35 | 显示全部楼层

回复 5楼 的帖子

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

[ 本帖最后由 dhp 于 2008-1-17 09:37 编辑 ]
发表于 2008-1-17 09:53 | 显示全部楼层


外循环已经给定了ξ的值,这个已经没有关系了啊
发表于 2008-1-17 09:58 | 显示全部楼层
原帖由 dhp 于 2008-1-17 09:35 发表
对,int是符号积分。不过,它也能算出数值,虽然有warning。我对int和quad的区别理解还不深。   
thx


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

quad 等是数值积分,根据simpson,gauss等数值求积的方法来计算,可以参考数值积分的书
 楼主| 发表于 2008-1-17 10:00 | 显示全部楼层

回复 7楼 的帖子

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

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

[ 本帖最后由 sigma665 于 2008-1-17 10:06 编辑 ]
 楼主| 发表于 2008-1-18 04:45 | 显示全部楼层

继续昨天的题目请教数值积分!

目前我完成了整个数值积分过程,与刚开始的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 编辑 ]
发表于 2008-1-18 10:14 | 显示全部楼层

回复 楼主 的帖子

singularity possible

积分表达式可能是奇异的...
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-18 11:41 , Processed in 0.065637 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表