声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4350|回复: 3

[编程技巧] 二分法matlab求解的程序 --给mm

[复制链接]
发表于 2005-5-27 20:39 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
发信人: zjliu (秋天的萝卜), 信区: Matlab
标 题: Re: 请教,这个方程怎么解
发信站: 哈工大紫丁香 (Fri May 27 20:01:33 2005), 转信
这是二分法求解的程序,结果x的值显示在figure上

  1. function lushan
  2. %
  3. %
  4. cos(alpha)=-3.8*w/sqrt((2.7*f)^2+(3.8w)^2);
  5. %
  6. cos(beta)=-2.7*f/sqrt((2.7*f)^2+(3.8w)^2);
  7. % f=38+3.8*x*cos(alpha);
  8. %
  9. w=54+2.7*x*cos(beta);
  10. % f*w-1140=0;
  11. close
  12. all
  13. figure;
  14. ezplot('2.7^2*f*(f-38)-3.8^2*w*(w-54)',[-1000,1000]);
  15. p=get(gca,'Children');
  16. set(p,'color','r');hold
  17. on;
  18. ezplot('w*f-1140',[-1000,1000]);
  19. ha=gca;figure;copyobj(ha,gcf)
  20. xlim([-60,40]);fl=10;fr=40;
  21. W=solve('w*(w-54)-a=0','w');
  22. while
  23. abs(fl-fr)>1e-4;
  24. f=(fl+fr)/2;
  25. w1=1140/f;
  26. a=2.7^2*f*(f-38)/3.8^2;
  27. w2=subs(W(1),a);
  28. if
  29. w1>w2;
  30. fl=f;
  31. else
  32. fr=f;
  33. end
  34. w=(w1+w2)/2;
  35. end
  36. text(f-25,w+300,
  37. ['(w=',num2str(w),',f=',num2str(f),')']);
  38. alphaa=-3.8*w/sqrt((2.7*f)^2+(3.8*w)^2);
  39. x=(f-38)
  40. /(3.8*cos(cos(alphaa)));
  41. text(f-25,w+200,['解一:
  42. x=',num2str(f)],'fontsize',16,...
  43. 'color','r');
  44. plot(f,w,'*');plot(f,w,'s');
  45. %%%%%%%%%%%
  46. 第二个解
  47. fl=-60;fr=-30;
  48. W=solve('w*(w-54)-a=0','w');< BR>while
  49. abs(fl-fr)>1e-4;
  50. f=(fl+fr)/2;
  51. w1=1140/f;
  52. a=2.7^2*f*(f-38)/3.8^2;
  53. w2=subs(W(2),a);
  54. if
  55. w1>w2;
  56. fl=f;
  57. else
  58. fr=f;
  59. end
  60. w=(w1+w2)/2;
  61. end
  62. text(-50,-200,
  63. ['(w=',num2str(w),',f=',num2str(f),')']);
  64. alphaa=-3.8*w/sqrt((2.7*f)^2+(3.8*w)^2);
  65. x=(f-38)
  66. /(3.8*cos(cos(alphaa)));
  67. text(-50,-400,['解二:
  68. x=',num2str(f)],'fontsize',16,...
  69. 'color','r');
  70. plot(f,w,'*');plot(f,w,'s');
复制代码


回复
分享到:

使用道具 举报

发表于 2005-5-31 17:46 | 显示全部楼层
心灯抢我的程序,kick
给哪个mm阿?
发表于 2005-5-31 18:33 | 显示全部楼层

上传两篇模态的文章供大家参考

哈,上面不已经注明你的信息了么?:)
你声明不能转载了么?帮你宣传一下,没让你请客就是好的了。
发表于 2011-11-12 12:22 | 显示全部楼层
1F的程序好像没转贴好
无权编辑, 有权的人帮忙编辑1F下, 再删此帖!
  1. function zza
  2. %
  3. % cos(alpha)=-3.8*w/sqrt((2.7*f)^2+(3.8*w)^2);
  4. % cos(beta)=-2.7*f/sqrt((2.7*f)^2+(3.8*w)^2);
  5. % f=38+3.8*x*cos(alpha);
  6. % w=54+2.7*x*cos(beta);
  7. % f*w-1140=0;
  8. close all; figure;
  9. ezplot('2.7^2*f*(f-38)-3.8^2*w*(w-54)',[-1000,1000]);
  10. p=get(gca,'Children'); set(p,'color','r'); hold on;
  11. ezplot('w*f-1140',[-1000,1000]);
  12. ha=gca; figure; copyobj(ha,gcf)
  13. xlim([-60,40]); fl=10; fr=40;
  14. W=solve('w*(w-54)-a=0','w');
  15. while abs(fl-fr)>1e-4
  16.    f=(fl+fr)/2; w1=1140/f; a=2.7^2*f*(f-38)/3.8^2; w2=subs(W(1),a);
  17.    if w1>w2, fl=f; else fr=f; end
  18.    w=(w1+w2)/2;
  19. end
  20. text(f-25,w+300, ['(w=',num2str(w),',f=',num2str(f),')']);
  21. alphaa=-3.8*w/sqrt((2.7*f)^2+(3.8*w)^2);
  22. x=(f-38)/(3.8*cos(cos(alphaa)));
  23. text(f-25,w+200,['解一: x=',num2str(f)],'fontsize',16,'color','r');
  24. plot(f,w,'*');plot(f,w,'s');
  25. %%%%%%%%%%% 第二个解
  26. fl=-60; fr=-30;
  27. W=solve('w*(w-54)-a=0','w');
  28. while abs(fl-fr)>1e-4,
  29.    f=(fl+fr)/2; w1=1140/f; a=2.7^2*f*(f-38)/3.8^2; w2=subs(W(2),a);
  30.    if w1>w2, fl=f; else fr=f; end
  31.    w=(w1+w2)/2;
  32. end
  33. text(-50,-200, ['(w=',num2str(w),',f=',num2str(f),')']);
  34. alphaa=-3.8*w/sqrt((2.7*f)^2+(3.8*w)^2);
  35. x=(f-38)/(3.8*cos(cos(alphaa)));
  36. text(-50,-400,['解二: x=',num2str(f)],'fontsize',16,'color','r');
  37. plot(f,w,'*');plot(f,w,'s');
复制代码

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-16 06:00 , Processed in 0.053924 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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