subwaylo 发表于 2015-11-17 13:52

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

本帖最后由 subwaylo 于 2015-11-17 13:56 编辑


function F = myfun1(x,a)
mu=x(1);
delte=x(2);
alpha=a;
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));
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));
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));
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)));

F=;
x0 = ;         % 取初值
options = optimset('Display','off');
a=0:0.01:5;             % 变量取值范围
for i=1:1:length(a)
kk=a(i);
x = fsolve(@myfun1, x0, options,kk);%求解非线性方程组
x1(i)=x(1); x2(i)=x(2);
end
plot(a,x1,'-b',a,x2,'-r'); xlabel('k'); legend('x1','x2')

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

弗朗索瓦 发表于 2015-11-18 10:10

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

subwaylo 发表于 2015-11-22 11:46

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

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

subwaylo 发表于 2015-11-22 13:27


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

弗朗索瓦 发表于 2015-11-26 10:08

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

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

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

弗朗索瓦 发表于 2015-11-26 10:13

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

大概的思路是
k = 0;
x0 = ;
for i = 1:2
    x = fsolve(@myfun1, x0, options,kk);
    x0 = x;
    k=0.001
end
页: [1]
查看完整版本: fsolve数值求解一直运行没有结果怎么回事啊