查看: 266|回复: 6|关注: 0

[已解决] (暂未解决)vpasolve求解范围太大则得不到结果

[复制链接]

新手

40 麦片

财富积分


050


10

主题

55

帖子

1

最佳答案
10 财富积分
本帖最后由 南海鱼游 于 2020-5-22 20:18 编辑

使用vpasolve求解如下方程的零点,但当求解范围设置太大确往往得不到结果,即如下代码中,求解范围为[0,1e3]可以得到正确结果,而若范围为[0,1e6],则返回结果为空:

  1. b=500;c=1;h=10;
  2. syms x y
  3. y(x)=(b-x)./c.^2.*h.*exp(-(x-b).^2./(2.*c.^2))-(1e-6);
  4. vpasolve(y,[0,1e6])
  5. vpasolve(y,[0,1e3])
复制代码
运行结果:
  1. ans =
  2. Empty sym: 0-by-1

  3. ans =
  4. 494.01547450426081642041960218714
复制代码
实际过程中往往并不确定解所在的范围,于此向各位寻求建议。

排除疑问1:
曾以为是所给范围的端点的函数值符号相同导致找不到解,但运行以下代码,发现并不是:代码中函数y(x)在 0,1e3,1e6 的值均为负值。
  1. b=500;c=1;h=10;
  2. syms x y
  3. y(x)=(b-x)./c.^2.*h.*exp(-(x-b).^2./(2.*c.^2))-(1e-6);
  4. vpa([y(0),y(1e5)])
  5. vpasolve(y,[0,1e5])
  6. vpa([y(0),y(1e3)])
  7. vpasolve(y,[0,1e3])
复制代码
运行结果:
  1. ans =
  2. [ -0.000001, -0.000001]

  3. ans =
  4. Empty sym: 0-by-1

  5. ans =
  6. [ -0.000001, -0.000001]

  7. ans =
  8. 494.01547450426081642041960218714
复制代码






回复主题 已获打赏: 0 积分

举报

新手

40 麦片

财富积分


050


10

主题

55

帖子

1

最佳答案
 楼主| 发表于 2020-5-21 21:25:03 | 显示全部楼层
本帖最后由 南海鱼游 于 2020-5-22 20:01 编辑

尝试1:

想办法从公式推导上得到一个合适的解的范围或解的临近点,但没有实现;只确定解大于0;
暂时的尝试解决办法是选择一个参考值h,然后在一定范围内产生一个随机数,两者相乘,则以(0,h*rand)为解的搜寻范围,若得不到,则循环,直到循环10次:

  1. b=500;c=1;h=10;
  2. syms x y
  3. y(x)=(b-x)./c.^2.*h.*exp(-(x-b).^2./(2.*c.^2))-(1e-6);
  4. d_v=vpasolve(y,[0,1e6]);d_v=1;
  5. deltaD=y(d_v);
  6. if isempty(d_v)
  7.     deltaD=1;
  8. end
  9. countLoop=0;
  10. while ~(abs(deltaD)<(1e-10) || countLoop==10)
  11.     ra=1e-6;rb=1e2;
  12.     r=ra+(rb-ra).*rand(1,1);
  13.     d_v=vpasolve(y,[0,h.*r]);
  14.     deltaD=y(d_v);
  15.     if isempty(d_v)
  16.         deltaD=1;
  17.     end
  18.     countLoop=countLoop+1
  19. end
  20. [deltaD;d_v]
复制代码

然而这种方法也并不总是能解决问题,有很大的不确定性。



回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

中级

1151 麦片

财富积分


5001500


0

主题

2707

帖子

248

最佳答案
  • 关注者: 164
发表于 2020-5-21 22:42:41 | 显示全部楼层
参考下下面的解:
x: 494.015474504261
回复此楼 已获打赏: 0 积分

举报

新手

40 麦片

财富积分


050


10

主题

55

帖子

1

最佳答案
 楼主| 发表于 2020-5-21 22:55:08 | 显示全部楼层
shihe 发表于 2020-5-21 22:42
参考下下面的解:
x: 494.015474504261

这个解通过修改缩小求解范围可以得到,如上所述,vpasolve(y,[0,1e6]),得不到结果,而 vpasolve(y,[0,1e3])  确可以得到结果;
所以我的问题是如何避免这种因为求解范围太大而找不到解的问题。
回复此楼 已获打赏: 0 积分

举报

新手

40 麦片

财富积分


050


10

主题

55

帖子

1

最佳答案
 楼主| 发表于 2020-5-22 20:08:40 | 显示全部楼层
尝试2:

关键在于不知道指定函数解的确定范围,尝试通过观察图像来手动确定求解范围;
以下代码在指定的范围求解失败后,则画出尽量大的范围的函数图像,然后通过观察得的解的范围输入弹出的对话框,而后根据此范围求解;再次求解失败后画出上次指定范围的函数图像和对话框。
这种方法好处是人工干预求解范围,只要有耐心,一定可以得到解;缺点是,效率太低。:Q

  1. b=500;c=1;h=10;
  2. syms x y
  3. y(x)=(b-x)./c.^2.*h.*exp(-(x-b).^2./(2.*c.^2))-(1e-6);
  4. x0=vpasolve(y,[0,1e6]);
  5. deltaD=y(x0);
  6. if isempty(x0)
  7.     deltaD=1;
  8. end
  9. countLoop=0;
  10. %%%
  11. %Try some other methods before the inputdlg
  12. %%%
  13. tryLeft=-1e3;tryRigth=1e3;
  14. while abs(deltaD)>(1e-10)
  15.     plotFuncOn(tryLeft,tryRigth,f1_d,1);
  16.     prompt = {'x Start:','x End:'};
  17.     title = 'Input Scale ';
  18.     dims = [1 35];
  19.     definput = {'-1000','1000'};
  20.     xCell = inputdlg(prompt,title,dims,definput);
  21.     tryLeft=str2double(xCell{1});tryRight=str2double(xCell{2});
  22.     x0=vpasolve(y,[tryLeft,tryRight]);
  23.     deltaD=vpa(y(x0));
  24.     if isempty(deltaD)
  25.         deltaD=1;
  26.     end
  27.     countLoop=countLoop+1;
  28.     fprintf('\ncountLoop=%d\tdeltaD=%f',countLoop,deltaD);
  29. end
  30. close(1);
  31. x0

  32. function [h_fig]=plotFuncOn(initX,endX,funcP,figN)
  33. var='X_';
  34. syms X_
  35. Y(X_)=subs(vpa(funcP),var);
  36. xPoints=linspace(initX,endX,800);
  37. % if ~isempty(find(xPoints==-250))
  38. %     pause
  39. % end
  40. yPoints=Y(xPoints);
  41. y1=min(real(yPoints));
  42. y2=max(real(yPoints));
  43. h_fig=figure(figN);
  44. % axis([initX,endX,double(y1),double(y2)]);
  45. plt=plot(xPoints,yPoints);
  46. plt.LineWidth=1.5;
  47. % plot(xPoints,yPoints,'.');
  48. grid on;
  49. end
复制代码


回复此楼 已获打赏: 6 积分

举报

论坛优秀回答者

权威

4341 麦片

财富积分



4

主题

4396

帖子

940

最佳答案
  • 关注者: 213
发表于 2020-5-27 09:34:11 | 显示全部楼层 |此回复为最佳答案
南海鱼游 发表于 2020-5-21 22:55
这个解通过修改缩小求解范围可以得到,如上所述,vpasolve(y,[0,1e6]),得不到结果,而 vpasolve(y,[0,1e ...

所以我的问题是如何避免这种因为求解范围太大而找不到解的问题。

你为何这样想?过大的搜索空间,一般得不到结果!这很好理解啊!
你在中国找个人与在某个镇找个人,时间与可能性一样吗?
特别对于比较特别的方程!
这道题,再约束一下还可以得到另一个根!
499.9999998999999999999995
回复此楼 已获打赏: 0 积分

举报

新手

40 麦片

财富积分


050


10

主题

55

帖子

1

最佳答案
 楼主| 发表于 2020-5-29 11:41:13 | 显示全部楼层
本帖最后由 南海鱼游 于 2020-5-29 11:42 编辑
maple1314168 发表于 2020-5-27 09:34
所以我的问题是如何避免这种因为求解范围太大而找不到解的问题。

你为何这样想?过大的搜索空间,一般得 ...

是的,求解范围大意味着 计算量大和计算时间长;
但从理论上讲,更大的范围比更小的范围得到解的可能性更大,不是吗?
这就牵涉到vpasolve()的算法问题,但源码没读太懂;
就自己写了一个算法,用于求解单调函数的零点(指定曲线连续的边界,不计算边界点);
代码如下;但使用时有时会遇到除零报错的问题

  1. clear;clc;
  2. syms x y
  3. y=1/(2.*x^7+x^3+2)+100;
  4. y2=x+3;
  5. [x0,dy]=solveMono(y,[-100,0])

  6. function [x0,y0]=solveMono(fHd,limVec)
  7. nGrid=1000;%Search Grid;
  8. p=0.5;%march Infactor
  9. err=1e-12;% zero precision
  10. dErr=1e-20;%dy precision
  11. x0=mean(limVec);
  12. syms X_
  13. y(X_)=subs(fHd,X_);
  14. yp(X_)=diff(y,X_);
  15. y0=y(x0);yp0=yp(x0);
  16. % y0 + yp0 + stepSign -
  17. % y0 - yp0 - stepSign -
  18. % y0 + yp0 - stepSign +
  19. % y0 - yp0 + stepSign +
  20. count=0;
  21. stepSign=logical(y0.*yp0<0)-logical(y0.*yp0>0);
  22. lim=limVec((stepSign-(-1))/2+1);
  23. dy=1+err;
  24. figure;
  25. while abs(y0)>err && abs(dy)>err
  26.     xVec1=linspace(x0,lim,nGrid);
  27.     xVec=xVec1(1:end-1);
  28.     yVec=y(xVec);
  29.     productY=yVec(1:end-1).*yVec(2:end);
  30.     negIndex=find(productY<0);
  31.     if isempty(negIndex)
  32.         x0Last=xVec(floor(nGrid.*p));y0Last=yVec(floor(nGrid.*p));
  33.         x0=xVec(floor(nGrid.*p)+1);y0=yVec(floor(nGrid.*p)+1);
  34.     else
  35.         x0Last=xVec(negIndex(1));y0Last=yVec(negIndex(1));
  36.         x0=xVec(negIndex(1)+1);y0=yVec(negIndex(1)+1);
  37.     end
  38.     count=count+1;
  39.     dy=y0-y0Last;
  40.     fprintf('\nsolveMono():%d,%.12f,%.12f',count,x0,y0);
  41. %     vpa([count,x0,y0])
  42.     if count==1
  43.         subplot(2,2,1);hold on;plot(x0Last,count,'o');xlabel('x_0');ylabel('count');drawnow;title('Search Contrary Sign');grid on;
  44.         subplot(2,2,3);hold on;plot(count,y0Last,'o');xlabel('count');ylabel('y_0');drawnow;grid on;
  45.     end
  46.     subplot(2,2,1);hold on;plot(x0,count,'o');xlabel('x_0');ylabel('count');drawnow;title('Search Contrary Sign');grid on;
  47.     subplot(2,2,3);hold on;plot(count,y0,'o');xlabel('count');ylabel('y_0');drawnow;grid on;
  48.     if y0*y0Last<0
  49.         if y0<y0Last
  50.             xLeft=x0;xRight=x0Last;
  51.             yLeft=y0;%yRight=y0Last;
  52.         elseif y0>=y0Last
  53.             xLeft=x0Last;xRight=x0;
  54.             yLeft=y0Last;%yRight=y0;
  55.         end
  56.         while abs(yLeft)>err && abs(dy)>dErr
  57.             xVec=linspace(xLeft,xRight,nGrid);
  58.             yVec=y(xVec);
  59.             productY=yVec(1:end-1).*yVec(2:end);
  60.             negIndex=find(productY<0);
  61.             yLeftLast=yLeft;
  62.             xLeft=xVec(negIndex(1));yLeft=yVec(negIndex(1));
  63.             xRight=xVec(negIndex(1)+1);%yRight=yVec(negIndex(1)+1);
  64.             dy=yLeft-yLeftLast;
  65.             count=count+1;
  66.             fprintf('\nsolveMono():%d,%.12f,%.12f',count,xLeft,yLeft);
  67. %             vpa([count,xLeft,yLeft])
  68.             subplot(2,2,2);hold on;plot(xLeft,count,'*');xlabel('x_0');ylabel('count');drawnow;title('Search Zero Point');grid on;
  69.             subplot(2,2,4);hold on;plot(count,yLeft,'*');xlabel('count');ylabel('y_0');drawnow;grid on;
  70.         end
  71.         x0=xLeft;y0=double(yLeft);
  72.     end
  73. end
  74. close(gcf)
  75. end
复制代码
回复此楼 已获打赏: 0 积分

举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

快速回复 返回顶部 返回列表