matlab对复杂积分函数作图问题

悠若谷

积分上限x0可以通过n（x0）代入上面第一个方程求出来
结果是一个关于delta-n跟D的式子
当积分项跟积分上限同时含有未知数D的时候
程序应该怎么来实现？
哪位有想法的给小妹提供个思路吧！多谢了！

悠若谷

下面是我的完成一半的程序，不知道下面该怎么进展了
function y=mytdraw(x)
syms y;
N0=1.554336; %赋值，x-deltan y-D
k=2*pi/(632.8e-9);
B0=k*N0;
ns=1.521145;
function f1=f1(z) %积分项用内嵌函数表示
f1=k*sqrt((ns+x.*erfc(z./y)).^2-N0^2);%这里已经用到y了，这么写可以吗
end
R1=quad(@f1,0,y*finverse_erfc((N0-ns)/x));%积分 x0=y*finverse_erfc((N0-ns)/x))里边有y不能这么用了但不知道咋改
eq1=R1-pi/4-atan(sqrt((B0^2-k^2)/(k^2*(ns+x)^2-B0^2))-x*(k^2)*(ns+x)/(pi^0.5*y*(k^2*(ns+x)^2-B0^2)^1.5));
y=fzero(eq1,[5e-7,5e-5]);  
end

comman window
clear all
x=linspace(0.01,0.08,801);%x赋值 挨个画图
y=zeros(size(x));
for i=1:length(x)
y(i)=mytdraw(x(i));
end
plot(x,y)

rocwoods

1. 
   
2. function y=mytdraw(x)
   
3. N0=1.554336; %赋值，x-deltan y-D
   
4. k=2*pi/(632.8e-9);
   
5. B0=k*N0;
   
6. ns=1.521145;
   
7. % function f1=f1(z) %积分项用内嵌函数表示
   
8. % f1=k*sqrt((ns+x.*erfc(z./y)).^2-N0^2);%这里已经用到y了，这么写可以吗
   
9. % end
   
10. eq1 = @(y) quadl(@(z) k*sqrt((ns+x.*erfc(z./y)).^2-N0^2),0,y*erfcinv((N0-ns)/x))-...
    
11.     pi/4-atan(sqrt((B0^2-k^2)/(k^2*(ns+x)^2-B0^2))-x*(k^2)*(ns+x)/(pi^0.5*y*(k^2*(ns+x)^2-B0^2)^1.5));
    
12. y=fzero(eq1,0.000001);  
    
13. end
    

复制代码

想来想去，你这个例子不错，于是给出完整代码解决方法，可以解决这样一类问题求解方法，方便需要的人参考。上面代码格式应该没有问题，需要说明的是，我计算finverse_erfc((N0-ns)/0.01)时候会报错，所以用MATLAB自带的erfcinv函数替代了finverse_erfc了。楼主需要的应该是erfcinv函数吧？
按以上代码在运行clear all
x=linspace(0.01,0.08,801);%x赋值 挨个画图
y=zeros(size(x));
for i=1:length(x)
y(i)=mytdraw(x(i));
end
plot(x,y)
也会报错，但是x=linspace(0.04,0.08,800)就不会报错。原因可能和参数x范围有关。楼主可以根据原方程物理意义，看看需要将x定为多少。
下图是我x按linspace(0.04,0.08,800)画出来的delta-D的图,横轴就是x(delta).
不懂你的方程物理意义，直觉上感觉delta的范围是不是有问题？

[ 本帖最后由 rocwoods 于 2009-6-22 01:18 编辑 ]

悠若谷

delta-n的范围是我贴错了 要大于0.0339的 前辈好强哦！
小声得说一句，我才知道matlab自带erfcinv函数。。。
我用萝卜前辈写的finverse_erfc，运行结果有些值能计算有些不能，昨天还想放弃这个程序呢
用erfcinv就没问题了！
太谢谢rocwoods前辈了！~