螺杆转子刀具设计（2） matlab 解法

shouce

螺杆转子刀具设计 数学模型
. U! M$ u# T  v; {9 @
/ f! ~" S2 X' T解如下方程：1679999.7*sin(x)-1679997.5*cos(y)*sin(x)-1499997.5*sin(y)*cos(x)-75000*sin(x)^3-179999.4000005*y*cos(y)*cos(x)+179999.4000005*y*sin(y)*sin(x)+75000*cos(y)*cos(x)*sin(x)=0
求当x=-1.361时  y=?  
>> syms x y
>> f=1679999.7*sin(x)-1679997.5*cos(y)*sin(x)-1499997.5*sin(y)*cos(x)-75000*sin(x)^3-179999.4000005*y*cos(y)*cos

6 P2 s( y. s% `! B1 x3 \5 v(x)+179999.4000005*y*sin(y)*sin(x)+75000*cos(y)*cos(x)*sin(x)
0 Y0 X  u$ s6 G, @+ t( V) ^2 V$ x
f =
; |+ i7 h1 ~& ~) }) l
(7215543768789811*sin(x))/4294967296 - (2999995*cos(x)*sin(y))/2 - (3359995*cos(y)*sin(x))/2 - 75000*sin(x)^3 + 75000*cos(x)*cos(y)

" M8 U3 {' U8 S* J* _9 y*sin(x) - (6184732290414159*y*cos(x)*cos(y))/34359738368 + (6184732290414159*y*sin(x)*sin(y))/34359738368

>>  subs(f, x, -1.361)

ans =
4 ]! s) r( \9 i0 C% ]+ W3 ^3 V
(3359995*sin(1361/1000)*cos(y))/2 - (2999995*cos(1361/1000)*sin(y))/2 - (7215543768789811*sin(1361/1000))/4294967296 + 75000*sin
4 _# g! {2 n4 F* ]- {
1 X: U6 Z+ u8 w(1361/1000)^3 - 75000*cos(1361/1000)*sin(1361/1000)*cos(y) - (6184732290414159*y*cos(1361/1000)*cos(y))/34359738368 -

(6184732290414159*y*sin(1361/1000)*sin(y))/34359738368


>>  a=0; b=1;
eps1=1e-8;eps2=1e-8;
N=300;
f=@ (x) ((3359995*sin(1361/1000)*cos(x))/2 - (2999995*cos(1361/1000)*sin(x))/2 - (7215543768789811*sin(1361/1000))/4294967296 + 75000*sin(1361/1000)^3- 75000*cos(1361/1000)*sin(1361/1000)*cos(x) - (6184732290414159*x*cos(1361/1000)*cos(x))/34359738368 -(6184732290414159*x*sin(1361/1000)*sin(x))/34359738368);
Hfun=@Bisection;
  \* R- ~" H! M" t[k,x,f_value]= feval(Hfun, f,a,b,eps1,eps2,N);
  t# u- R! @' M- y# p$ A
运行结果
2 a$ B  d1 Q# y0 d! j, Nk ,a ,b ,x, f
( t* }, c- \& [* z+ Y& Y  1, 0.000000000, 1.000000000, 0.500000000,-352805.622314164,
) u+ g  o' J7 g2 `4 s  2, 0.000000000, 0.500000000, 0.250000000,-92968.331400711,
  3, 0.000000000, 0.250000000, 0.125000000,-4146.780462183,
% c7 c% `7 U" ~3 F8 F0 n/ J  4, 0.000000000, 0.125000000, 0.062500000,29178.936016433,
  5, 0.062500000, 0.125000000, 0.093750000,13458.364371693,
  6, 0.093750000, 0.125000000, 0.109375000,4890.124291139,
3 Q# z+ a8 E; o  7, 0.109375000, 0.125000000, 0.117187500,430.092693351,
9 w- L0 x* r. M, S1 D9 S! t2 J* B  8, 0.117187500, 0.125000000, 0.121093750,-1843.759437944,
2 X/ G: R6 Y: ^5 m8 L5 P2 i  9, 0.117187500, 0.121093750, 0.119140625,-703.184646645,
10, 0.117187500, 0.119140625, 0.118164063,-135.633470273,
11, 0.117187500, 0.118164063, 0.117675781,147.457778641,
12, 0.117675781, 0.118164063, 0.117919922,5.969190901,
: A  R% u. M" s+ A 13, 0.117919922, 0.118164063, 0.118041992,-64.817881139,
4 T$ \& _/ s. z, p# o% _/ Y 14, 0.117919922, 0.118041992, 0.117980957,-29.420780403,
15, 0.117919922, 0.117980957, 0.117950439,-11.724903562,
16, 0.117919922, 0.117950439, 0.117935181,-2.877633532,
9 y9 o# n# ^5 {2 R( T 17, 0.117919922, 0.117935181, 0.117927551,1.545834384,
18, 0.117927551, 0.117935181, 0.117931366,-0.665885649,
19, 0.117927551, 0.117931366, 0.117929459,0.439977849,
20, 0.117929459, 0.117931366, 0.117930412,-0.112953030,
  A$ y& s# Y$ |! G% w0 c2 s 21, 0.117929459, 0.117930412, 0.117929935,0.163512627,
- a6 K  }% Q- K 22, 0.117929935, 0.117930412, 0.117930174,0.025279853,
4 ~- K2 {" l' s6 h 23, 0.117930174, 0.117930412, 0.117930293,-0.043836575,
24, 0.117930174, 0.117930293, 0.117930233,-0.009278357,
& W$ T- B. u+ P6 z 25, 0.117930174, 0.117930233, 0.117930204,0.008000749,
/ s! O& H# D3 n: p0 B! j 26, 0.117930204, 0.117930233, 0.117930219,-0.000638804,
/ d3 G, L7 R* O) O; R" r* b+ T 27, 0.117930204, 0.117930219, 0.117930211,0.003680972,
8 F2 _- G3 W3 @4 [. K; k% k5 I>>
; W. Y/ x; c; X# m& I7 n  w* p( _
2 i& y9 \9 O& U" O! Q; S9 [结论 当x=-1.361时  y=0.117930219
按此方法   可依次  得出150对x与y的值
6 \' B! E! @& r: B3 F

单车居士

西交大学生？？问邢老师

明月山河

你这个是代数方程？

cosxuan

算法？

shouce

cosxuan 发表于 2015-6-4 19:45
算法？

二分法        
还有一个程序
% Bisection.m
function[k,x,f_value]=Bisection(f,a,b,eps1,eps2,N)
fprintf('k ,a ,b ,x, f\n')
3 L' ~* u; @, ^- U7 _8 e7 y& hfor k=1:N
    x=(a+b)/2;
    f_value=f(x);
& x5 B. I4 `3 n( w8 R    fprintf( '%3d, %10.9f, %10.9f, %10.9f,%10.9f,\n'...
1 U6 A) a& n: r      ,k ,a ,b ,x, f_value)
    if abs(f_value)< eps1||0.5*(b-a)<eps2
    return
' M8 P8 S2 A- d5 o, k2 Z8 [    else
        if f(x)*f(a)<0
          b=x;
; V/ C4 @) k1 s( D. h        else
& Z4 R) R8 |2 b) k          a=x;
0 z$ z+ H+ R& d: X        end
+ B9 y6 e3 f; Z. b        if k== N
6 I/ A4 I0 G! W  H          warning ('算法超出最大迭代数!')   
        end
    end
- S) Y; L* N8 ]; x7 X% l4 R: ?+ O5 P

shouce

单车居士 发表于 2015-6-4 12:41
  W! H9 J9 c: p西交大学生？？问邢老师

7 z" E1 Q" L5 {: \- _9 G这点问题   自己就可以解决
5 Q  `) l9 t. f  |. W. x1 g8 ]

shouce

明月山河 发表于 2015-6-4 19:12
# F" q9 ]( t, }# S0 F( y# R你这个是代数方程？

是的呀      不过matlab写出来的方程  括号很多  不好看   
4 L+ Z! @& j; G! J" [0 o

shouce

明月山河 发表于 2015-6-4 19:12
( S) J$ L; O! B; W你这个是代数方程？

% S: j- B5 f  p' V- n' [是超越方程      
6 U/ n6 T3 B% G
9 _/ d6 p' K! R3 q7 m# K# }: w  ]
[chāo yuè fāng chéng
1 y6 g* z: A- `7 }# l3 t