悠若谷 发表于 2009-6-9 19:29

这个积分方程matlab应该怎么求?


xco与xc1是两个未知数的函数
方程左边的积分部分尝试用quadl函数来表达
fun1=inline('k.*sqrt((1.521145+x(1).*(1-erf(x/(2.*(sqrt(x(2).*t))))))^2-N1^2)');
R1=quadl(fun1,0,x0,1e-8);
fun2=inline('k.*sqrt((1.521145+x(1).*(1-erf(x/(2.*(sqrt(x(2).*t))))))^2-N2^2)');
R2=quadl(fun2,0,x1,1e-8);

但是运行出错
??? Error using ==> inline.feval
Not enough inputs to inline function.
Error in ==> quadl at 64
y = feval(f,x,varargin{:}); y = y(:).';
Error in ==> only at 14
R1=quadl(fun1,0,x0,1e-8);
Error in ==> fsolve at 180
      fuser = feval(funfcn{3},x,varargin{:});

哪位能帮忙看看是怎么回事吗?

去掉inline函数 改用下面的程序
R1=quadl('k*sqrt((1.521145+x(1)*(1-erf(x/(2*(sqrt(x(2)*t))))))^2-N1^2)',0,x0,1e-8);
R2=quadl('k*sqrt((1.521145+x(1)*(1-erf(x/(2*(sqrt(x(2)*t))))))^2-N2^2)',0,x1,1e-8);
运行结果一样

[ 本帖最后由 ChaChing 于 2009-6-21 15:48 编辑 ]

ChaChing 发表于 2009-6-9 20:51

x0给什麽?

悠若谷 发表于 2009-6-10 13:42

x0在主函数fsolve函数里边给了初值的
'k*sqrt((1.521145+x(1)*(1-erf(x/(2*(sqrt(x(2)*t))))))^2-N1^2)'这个式子里边含有erf()函数
程序调试到这一步就停了,是不是对erf函数没法直接积分呢?那要怎么改

******only.m函数******
function y=only(x)
t=1200; N0=1.554336; N1=1.536857; ns=1.521145;
k=2*pi/(632.8e-9); B0=k*N0; B1=k*N1;
x0=2*finverse_erfc((N0-ns)/x(1))*sqrt(x(2)*t);
x1=2*finverse_erfc((N1-ns)/x(1))*sqrt(x(2)*t);
ierfcx1=int('k*sqrt((1.521145+x(1)*erfc(z/(2*(sqrt(x(2)*t)))))^2-N0^2)','z',0,x0);
ierfcx2=int('k*sqrt((1.521145+x(1)*erfc(z/(2*(sqrt(x(2)*t)))))^2-N1^2)','z',0,x1);
R1=subs(ierfcx1); R2=subs(ierfcx2);
tg1=sqrt((B0^2-k^2)/(k^2*(ns+x(1))^2-B0^2))-x(1)*(k^2)*(ns+x(1))/(2*sqrt(pi*x(2)*t)*(k^2*(ns+x(1))^2-B0^2)^1.5);
tg2=sqrt((B1^2-k^2)/(k^2*(ns+x(1))^2-B1^2))-x(1)*(k^2)*(ns+x(1))/(2*sqrt(pi*x(2)*t)*(k^2*(ns+x(1))^2-B1^2)^1.5);
eq1=R1-pi/4-atan(tg1); eq2=R2-5*pi/4-atan(tg2);
y=;
**********finverse_erfc.m函数********
function c=finverse_erfc(z);
% count inverse function of erfc
% author email of the program:% zjliu2001@163.com
a=600; b=0;
if z>1 | z<0;
   error([['    z takes a error value!'],...
         ['; please check it''s value!']]);
end
   
c=(a+b)/2;
while abs(erfc(c)-z)>1e-6;
   c=(a+b)/2;
   if erfc(c)-z<0, a=c; else b=c; end
end
********主函数******
x0=; options=optimset('display','iter');
=fsolve('only',x0,options)

[ 本帖最后由 ChaChing 于 2009-6-21 15:54 编辑 ]

rocwoods 发表于 2009-6-10 14:46

楼主还用了萝卜写的程序了啊?呵呵,你写的程序得需要较大的改动,晚上有时间可以帮你想想这个问题。

悠若谷 发表于 2009-6-10 15:58

楼上说的是finverse_erfc.m函数啊?
呵呵,在网上看到的就下下来了,感觉还挺好的!
嗯,谢谢recwoods~

调试时发现,运行到这一步的时候
ierfcx1=int('k*sqrt((1.521145+x(1)*erfc(z/(2*(sqrt(x(2)*t)))))^2-N0^2)','z',0,x0);
ierfcx2=int('k*sqrt((1.521145+x(1)*erfc(z/(2*(sqrt(x(2)*t)))))^2-N1^2)','z',0,x1);
fsolve设定的x初值没有传递进int积分函数
这个是为什么?应该怎么改呢?

[ 本帖最后由 ChaChing 于 2009-6-21 15:56 编辑 ]

rocwoods 发表于 2009-6-11 00:51

建议楼主only函数用nested function形式写,这样可以解决参数传递问题。楼主可以查看帮助文档或者google、baidu了解nested function用法。
关于向被积函数传参数,楼主可以参考这个帖子。
http://forum.vibunion.com/forum/thread-42369-1-1.html

悠若谷 发表于 2009-6-11 13:48

谢谢rocwoods~~
研究了下您说的nested function,
有个问题想问一下
我这里边需要编两个嵌套函数,s1跟s2是同级的。
那我可以这样使用吗?
function r = fatherfunction()
blah
   function s1 = childfunction()
    blah
   function s2 = childfunction()
   end
   end
end

悠若谷 发表于 2009-6-11 15:18

function f1=f1(z)
f1=sqrt(k^2.*(ns+x(1).*erfc(z/(2*(x(2)*t)^0.5))).^2-B0^2);
end
function f2=f2(z)
f2=sqrt(k^2.*(ns+x(1).*erfc(z/(2*(x(2)*t)^0.5))).^2-B1^2);
end
ierfcx1=quadl(@f1,0,x0); ierfcx2=quadl(@f2,0,x1);

程序运行到f1=sqrt(k^2.*(ns+x(1).*erfc(z/(2*(x(2)*t)^0.5))).^2-B0^2);就停住了

x0=; options=optimset('display','iter');
=fsolve('only',x0,options)

                                       Norm of      First-order   Trust-region
IterationFunc-count   f(x)          step         optimality    radius
   0          3         2789.92                     2.49e+013               1
MATLAB一直显示busy。。。

[ 本帖最后由 ChaChing 于 2009-6-21 16:00 编辑 ]

rocwoods 发表于 2009-6-11 22:46

回复 7楼 悠若谷 的帖子

你的blah是什么?在另一个帖子里如何求erfc(x)函数的平方的积分我已经给你回了,你可以参考下。你可以把你卡住的问题以最简单的例子呈现出来,不是光学专业的,不懂你那个方程各个物理量含义,方程看起来就费劲,不过我觉得你卡住的地方是可以用简单的例子说明白的。还有,工作后很少用QQ了,有问题可以在论坛上提,这样讨论的结果可以方便记录下来,同时方便以后可能遇到类似问题的人,这也是论坛的初衷之一。

[ 本帖最后由 ChaChing 于 2009-6-21 16:03 编辑 ]

悠若谷 发表于 2009-6-12 14:14

改了一下现在不是卡住了,设置断点单步运行的时候fsolve设定的初值传进去没有问题,第一,二次搜索很正常,但是到第三次搜索的时候程序就中断了。。。我了解的matlab知识不够用的,不知道是不是方程太复杂了不适合用fsolve。

                                       Norm of      First-order   Trust-region
IterationFunc-count   f(x)          step         optimality    radius
   0          3         9.66525                     2.03e+012               1
??? Input must be real.

Error in ==> erfc at 18
y = erfcore(x,1);

Error in ==> only>f1 at 14
f1=sqrt(k^2.*(ns+x(1)*erfc(z/(2*(x(2)*t)^0.5))).^2-B0^2);

Error in ==> quad at 62
y = f(x, varargin{:});

Error in ==> only at 16
R1=quad(@f1,0,x0);

Error in ==> optim\private\trustnleqn at 210
    F = feval(funfcn{3},x,varargin{:});

Error in ==> fsolve at 295
    =...

Error in ==> fmain at 3
=fsolve('only',x0,options)

[ 本帖最后由 悠若谷 于 2009-6-12 14:29 编辑 ]

悠若谷 发表于 2009-6-12 14:18

这个程序,如果把复杂的地方拿出来用简单的例子表示基本上运行是没问题的,但是都凑到一起就不行了。您要是有空的话帮我看看吧,这个程序的思想很简单,就是变量很多。
%******fmain函数*****
x0=;
options=optimset('display','iter');
=fsolve('only',x0,options)

%*****only.m********
function y=only(x)
t=1200;
N0=1.554336;
N1=1.536857;
ns=1.521145;
k=2*pi/(632.8e-9);
B0=k*N0;
B1=k*N1;
x0=2*finverse_erfc((N0-ns)/x(1))*sqrt(x(2)*t);
x1=2*finverse_erfc((N1-ns)/x(1))*sqrt(x(2)*t);
tg1=sqrt((B0^2-k^2)/(k^2*(ns+x(1))^2-B0^2))-x(1)*(k^2)*(ns+x(1))/(2*sqrt(pi*x(2)*t)*(k^2*(ns+x(1))^2-B0^2)^1.5);
tg2=sqrt((B1^2-k^2)/(k^2*(ns+x(1))^2-B1^2))-x(1)*(k^2)*(ns+x(1))/(2*sqrt(pi*x(2)*t)*(k^2*(ns+x(1))^2-B1^2)^1.5);
function f1=f1(z)
f1=sqrt(k^2.*(ns+x(1)*erfc(z/(2*(x(2)*t)^0.5))).^2-B0^2);
end
R1=quad(@f1,0,x0);
function f2=f2(z)
f2=sqrt(k^2.*(ns+x(1)*erfc(z/(2*(x(2)*t)^0.5))).^2-B1^2);
end
R2=quad(@f2,0,x1);
eq1=R1-pi/4-atan(tg1);
eq2=R2-5*pi/4-atan(tg2);
y=;
end

%***finverse_erfc.m******
function c=finverse_erfc(z);
% count inverse function of erfc
% author email of the program:% zjliu2001@163.com
a=600;
b=0;
if z>1 | z<0;
   error([['    z takes a error value!'],...
         ['; please check it''s value!']]);
end
   
c=(a+b)/2;
while abs(erfc(c)-z)>1e-6;
   c=(a+b)/2;
   if erfc(c)-z<0;
       a=c;
   else
       b=c;
   end
end
页: [1]
查看完整版本: 这个积分方程matlab应该怎么求?