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

shouce

螺杆转子刀具设计 数学模型
8 C, O' z  c( b' Q  L
解如下方程：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
# {5 y. J1 h. j  [# I求当x=-1.361时  y=?  
8 @9 h/ G) B; T( P. X" L; u>> syms x y
  [5 M9 m! _( v% R>> 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

$ d2 }6 F8 z- L1 _" v# u& b(x)+179999.4000005*y*sin(y)*sin(x)+75000*cos(y)*cos(x)*sin(x)
1 D. _; U8 `" F, z1 d
f =

% i' U; M. ]( j# [1 o% x(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)
* I5 X% f) w% G4 a% T7 u
*sin(x) - (6184732290414159*y*cos(x)*cos(y))/34359738368 + (6184732290414159*y*sin(x)*sin(y))/34359738368
2 h+ \/ r% i2 ~* R7 I0 u
>>  subs(f, x, -1.361)

ans =

: c' C4 t. E' E(3359995*sin(1361/1000)*cos(y))/2 - (2999995*cos(1361/1000)*sin(y))/2 - (7215543768789811*sin(1361/1000))/4294967296 + 75000*sin

(1361/1000)^3 - 75000*cos(1361/1000)*sin(1361/1000)*cos(y) - (6184732290414159*y*cos(1361/1000)*cos(y))/34359738368 -
' H5 v- p  j+ J
3 Q$ K; I" d7 l8 W4 F0 R(6184732290414159*y*sin(1361/1000)*sin(y))/34359738368


>>  a=0; b=1;
eps1=1e-8;eps2=1e-8;
) U$ N+ T! z2 Y- R7 F N=300;
4 I7 l. P, n/ ^  ~" Of=@ (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);
+ m9 B9 ^( F; r! x- @Hfun=@Bisection;
$ p. ~% m; h* ~- @, C[k,x,f_value]= feval(Hfun, f,a,b,eps1,eps2,N);

$ x5 a; l) x7 a1 `运行结果
5 L! a2 k. }- m" X$ Fk ,a ,b ,x, f
7 P6 u, n8 s$ t/ m- S  1, 0.000000000, 1.000000000, 0.500000000,-352805.622314164,
  2, 0.000000000, 0.500000000, 0.250000000,-92968.331400711,
  3, 0.000000000, 0.250000000, 0.125000000,-4146.780462183,
  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,
0 J0 _5 L2 A8 p5 T  @  7, 0.109375000, 0.125000000, 0.117187500,430.092693351,
  8, 0.117187500, 0.125000000, 0.121093750,-1843.759437944,
  9, 0.117187500, 0.121093750, 0.119140625,-703.184646645,
10, 0.117187500, 0.119140625, 0.118164063,-135.633470273,
. g( ]# I" _- h 11, 0.117187500, 0.118164063, 0.117675781,147.457778641,
12, 0.117675781, 0.118164063, 0.117919922,5.969190901,
* R. Q" n6 q) r+ {& C 13, 0.117919922, 0.118164063, 0.118041992,-64.817881139,
  w+ y$ Z! R9 E. _0 A: V3 j 14, 0.117919922, 0.118041992, 0.117980957,-29.420780403,
+ m7 D" S! P& ~% o$ O. R 15, 0.117919922, 0.117980957, 0.117950439,-11.724903562,
16, 0.117919922, 0.117950439, 0.117935181,-2.877633532,
17, 0.117919922, 0.117935181, 0.117927551,1.545834384,
18, 0.117927551, 0.117935181, 0.117931366,-0.665885649,
+ v* a# ]/ h+ o 19, 0.117927551, 0.117931366, 0.117929459,0.439977849,
5 u" A1 \2 a" M5 B# C 20, 0.117929459, 0.117931366, 0.117930412,-0.112953030,
! j) _) [4 n3 a$ K  T, r3 @ 21, 0.117929459, 0.117930412, 0.117929935,0.163512627,
- X! Z' b, W3 b9 ~  r# H 22, 0.117929935, 0.117930412, 0.117930174,0.025279853,
- O2 q0 |* s7 L. e; e/ ^7 t 23, 0.117930174, 0.117930412, 0.117930293,-0.043836575,
1 Q( |/ i  F" U6 { 24, 0.117930174, 0.117930293, 0.117930233,-0.009278357,
25, 0.117930174, 0.117930233, 0.117930204,0.008000749,
9 A; V6 L5 S7 d0 k# i 26, 0.117930204, 0.117930233, 0.117930219,-0.000638804,
27, 0.117930204, 0.117930219, 0.117930211,0.003680972,
0 Q1 s& k; o7 a4 L* `8 B$ {% y>>

" \- v7 a) Z$ C: |! P/ G结论 当x=-1.361时  y=0.117930219
9 m; V; q5 N) v, x& K  I  T按此方法   可依次  得出150对x与y的值
7 K/ Y6 D  K% C8 K) B

单车居士

西交大学生？？问邢老师

明月山河

你这个是代数方程？

cosxuan

算法？

shouce

cosxuan 发表于 2015-6-4 19:45
; M% e5 R- \& n" a! d+ z算法？

二分法        
. j. w: Z7 S- s5 u& i0 z% p还有一个程序
0 f6 H- q, E1 g4 c) W% Bisection.m
2 ^; B' ?: W/ d% s# j1 w' n7 Y- x' w, ifunction[k,x,f_value]=Bisection(f,a,b,eps1,eps2,N)
2 D" `/ ~- A8 r; ?9 afprintf('k ,a ,b ,x, f\n')
& Z5 [( Q' ^2 i% E3 mfor k=1:N
    x=(a+b)/2;
    f_value=f(x);
    fprintf( '%3d, %10.9f, %10.9f, %10.9f,%10.9f,\n'...
      ,k ,a ,b ,x, f_value)
    if abs(f_value)< eps1||0.5*(b-a)<eps2
    return
" a2 h' X' o! D3 e. p0 q" z1 O    else
        if f(x)*f(a)<0
2 p* b" }' r- Y3 F2 W          b=x;
        else
          a=x;
        end
$ H  u% Y8 M* u* z1 A! N  p        if k== N
/ F5 M) s; s" k! [9 j- }          warning ('算法超出最大迭代数!')   
        end
+ Z+ k* a2 H% _3 b    end

shouce

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

" y5 `& ]3 q! Q# f. E这点问题   自己就可以解决

shouce

明月山河 发表于 2015-6-4 19:12
你这个是代数方程？

是的呀      不过matlab写出来的方程  括号很多  不好看   
, x- L" I' r% ^0 L. @1 |
) K- ~7 s/ S. t! R: p

shouce

明月山河 发表于 2015-6-4 19:12
9 `: D% O) O6 K4 q% a1 @你这个是代数方程？

% V9 M; M1 @9 u" Q3 ~7 D6 t
) U' O" v, c. E) k& c( n是超越方程      
( P) T9 ~9 Q9 p3 J0 v
* ?0 ~& k( N  @; o5 o% Y7 z2 B
" U7 _6 s; g# X7 h' P: {0 Y7 v[chāo yuè fāng chéng
% a1 j9 M" G9 m/ u# h9 |