fsolve数值求解一直运行没有结果怎么回事啊

subwaylo

1. 
   

复制代码

1. function F = myfun1(x,a)
   
2. mu=x(1);
   
3. delte=x(2);
   
4. alpha=a;
   
5. i1z =- 2*ellipticE((-((alpha - mu)/(2*delte) - ((alpha - mu)^2/delte^2 + 1)^(1/2)/2)/((alpha - mu)^2/delte^2 + 1)^(1/2))^(1/2))*((alpha - mu)^2/delte^2 + 1)^(1/4) - (2*ellipticK((-((alpha - mu)/(2*delte) - ((alpha - mu)^2/delte^2 + 1)^(1/2)/2)/((alpha - mu)^2/delte^2 + 1)^(1/2))^(1/2)))/(((2*(alpha - mu))/delte - 2*((alpha - mu)^2/delte^2 + 1)^(1/2))*((alpha - mu)^2/delte^2 + 1)^(1/4)) - (ellipticK((-((alpha - mu)/(2*delte) - ((alpha - mu)^2/delte^2 + 1)^(1/2)/2)/((alpha - mu)^2/delte^2 + 1)^(1/2))^(1/2))*(alpha - mu))/(delte*((alpha - mu)^2/delte^2 + 1)^(1/4));
   
6. i1z1 =(2*ellipticK((((alpha + mu)/(2*delte) + ((alpha + mu)^2/delte^2 + 1)^(1/2)/2)/((alpha + mu)^2/delte^2 + 1)^(1/2))^(1/2)))/(((alpha + mu)^2/delte^2 + 1)^(1/4)*((2*(alpha + mu))/delte + 2*((alpha + mu)^2/delte^2 + 1)^(1/2))) - 2*ellipticE((((alpha + mu)/(2*delte) + ((alpha + mu)^2/delte^2 + 1)^(1/2)/2)/((alpha + mu)^2/delte^2 + 1)^(1/2))^(1/2))*((alpha + mu)^2/delte^2 + 1)^(1/4) + (ellipticK((((alpha + mu)/(2*delte) + ((alpha + mu)^2/delte^2 + 1)^(1/2)/2)/((alpha + mu)^2/delte^2 + 1)^(1/2))^(1/2))*(alpha + mu))/(delte*((alpha + mu)^2/delte^2 + 1)^(1/4));
   
7. i2z =ellipticK((-((alpha - mu)/(2*delte) - ((alpha - mu)^2/delte^2 + 1)^(1/2)/2)/((alpha - mu)^2/delte^2 + 1)^(1/2))^(1/2))/(3*((alpha - mu)^2/delte^2 + 1)^(1/4)) - (2*ellipticE((-((alpha - mu)/(2*delte) - ((alpha - mu)^2/delte^2 + 1)^(1/2)/2)/((alpha - mu)^2/delte^2 + 1)^(1/2))^(1/2))*((alpha - mu)^2/delte^2 + 1)^(1/4)*(alpha - mu))/(3*delte) - (2*ellipticK((-((alpha - mu)/(2*delte) - ((alpha - mu)^2/delte^2 + 1)^(1/2)/2)/((alpha - mu)^2/delte^2 + 1)^(1/2))^(1/2))*(alpha - mu))/(3*delte*((2*(alpha - mu))/delte - 2*((alpha - mu)^2/delte^2 + 1)^(1/2))*((alpha - mu)^2/delte^2 + 1)^(1/4));
   
8. i2z1 =ellipticK((((alpha + mu)/(2*delte) + ((alpha + mu)^2/delte^2 + 1)^(1/2)/2)/((alpha + mu)^2/delte^2 + 1)^(1/2))^(1/2))/(3*((alpha + mu)^2/delte^2 + 1)^(1/4)) + (2*ellipticE((((alpha + mu)/(2*delte) + ((alpha + mu)^2/delte^2 + 1)^(1/2)/2)/((alpha + mu)^2/delte^2 + 1)^(1/2))^(1/2))*((alpha + mu)^2/delte^2 + 1)^(1/4)*(alpha + mu))/(3*delte) - (2*ellipticK((((alpha + mu)/(2*delte) + ((alpha + mu)^2/delte^2 + 1)^(1/2)/2)/((alpha + mu)^2/delte^2 + 1)^(1/2))^(1/2))*(alpha + mu))/(3*delte*((alpha + mu)^2/delte^2 + 1)^(1/4)*((2*(alpha + mu))/delte + 2*((alpha + mu)^2/delte^2 + 1)^(1/2)));
   
9. 
   
10. F=[i1z+i1z1; i2z+i2z1-(4/3)*(delte)^(-3/2)];
    
11. x0 = [0.01; 0.01];           % 取初值
    
12. options = optimset('Display','off');
    
13. a=0:0.01:5;             % 变量取值范围[0 5]
    
14. for i=1:1:length(a)
    
15. kk=a(i);
    
16. x = fsolve(@myfun1, x0, options,kk);%求解非线性方程组
    
17. x1(i)=x(1); x2(i)=x(2);
    
18. end
    
19. plot(a,x1,'-b',a,x2,'-r'); xlabel('k'); legend('x1','x2')

复制代码

matlab新手，只是照着论坛上的程序模仿的，运行后直接卡死，是不是计算量太大，有没有好的解决办法，谢谢大家啦

弗朗索瓦

有可能是因为方程太复杂，一般情况下fsolve用于求解简单的非线性方程是没有问题的
如果方程比较特殊，建议根据方程的特征寻找适合的求解方法

subwaylo

弗朗索瓦 发表于 2015-11-18 10:10
有可能是因为方程太复杂，一般情况下fsolve用于求解简单的非线性方程是没有问题的
如果方程比较特殊，建议 ...

实在找不出其它方法啊，画是画出来了，但是初值影响很大，而且同一个初值我不知道能不能试用所有变量，请问还有什么好的办法能画出这个的

subwaylo

我想用第一次计算k=0的结果给第二次计算k=0.001赋初值该怎么循环编程啊，新手求教

弗朗索瓦

subwaylo 发表于 2015-11-22 11:46
实在找不出其它方法啊，画是画出来了，但是初值影响很大，而且同一个初值我不知道能不能试用所有变量，请 ...

fsolve属于局部最优算法，其结果敏感于初值是必然的，特别是对于比较方程
你可以尝试修改option中的选项看看是否能够改进

另外可以考虑使用1stOpt，看看求解的效果

弗朗索瓦

subwaylo 发表于 2015-11-22 13:27
我想用第一次计算k=0的结果给第二次计算k=0.001赋初值该怎么循环编程啊，新手求教

大概的思路是

1. k = 0;
   
2. x0 = [0.01; 0.01];
   
3. for i = 1:2
   
4.     x = fsolve(@myfun1, x0, options,kk);
   
5.     x0 = x;
   
6.     k=0.001
   
7. end

复制代码