这个积分方程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

x0给什麽?

悠若谷

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=[eq1;eq2]；
**********finverse_erfc.m函数********
function c=finverse_erfc(z);
% count inverse function of erfc[complementary error function]
% 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=[0.035;1e-12]; options=optimset('display','iter');
[x,fval]=fsolve('only',x0,options)

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

rocwoods

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

悠若谷

楼上说的是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

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

悠若谷

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

悠若谷

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=[0.035;1e-12]; options=optimset('display','iter');
[x,fval]=fsolve('only',x0,options)

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

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

rocwoods

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

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

悠若谷

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

                                         Norm of      First-order   Trust-region
Iteration  Func-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
    [x,FVAL,JACOB,EXITFLAG,OUTPUT,msg]=...

Error in ==> fmain at 3
[x,fval]=fsolve('only',x0,options)

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

悠若谷

这个程序，如果把复杂的地方拿出来用简单的例子表示基本上运行是没问题的，但是都凑到一起就不行了。您要是有空的话帮我看看吧，这个程序的思想很简单，就是变量很多。
%******fmain函数*****
x0=[0.0436;1e-14];
options=optimset('display','iter');
[x,fval]=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=[eq1;eq2];
end

%***finverse_erfc.m******
function c=finverse_erfc(z);
% count inverse function of erfc[complementary error function]
% 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