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

shouce

螺杆转子刀具设计 数学模型
5 x  Z( H2 t# X
6 U1 Q4 `+ x9 S; ^4 V解如下方程：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
0 r/ l  e+ l$ @" W# w>> 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
  `5 r; d* g) G# n$ a6 M
) d* l: p+ h) l- V6 i( }$ Y(x)+179999.4000005*y*sin(y)*sin(x)+75000*cos(y)*cos(x)*sin(x)

f =
: p5 d9 ]( h  t
3 O0 g' N$ ]' S) x) J4 J4 V. H(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)
0 r( j8 Q% l8 p+ {- A
6 Q- {: y' N# V8 D*sin(x) - (6184732290414159*y*cos(x)*cos(y))/34359738368 + (6184732290414159*y*sin(x)*sin(y))/34359738368
# a3 p! j$ f0 M: b4 w6 |
0 M0 l( E2 M  G( B>>  subs(f, x, -1.361)

ans =

! |. L5 n2 w- j5 R(3359995*sin(1361/1000)*cos(y))/2 - (2999995*cos(1361/1000)*sin(y))/2 - (7215543768789811*sin(1361/1000))/4294967296 + 75000*sin
. \$ e) E% Q0 A9 i8 E2 U7 i
(1361/1000)^3 - 75000*cos(1361/1000)*sin(1361/1000)*cos(y) - (6184732290414159*y*cos(1361/1000)*cos(y))/34359738368 -
3 A/ }; b' z# ]+ D$ T$ D
(6184732290414159*y*sin(1361/1000)*sin(y))/34359738368
0 k5 A3 \6 ]9 y, E
- ?% @* |( g3 H1 u% n, p3 r' t% w! D. D5 B
+ V0 J1 f, U, f>>  a=0; b=1;
eps1=1e-8;eps2=1e-8;
5 c: Y" H, h/ Q N=300;
# C+ A3 v# s) Tf=@ (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;
[k,x,f_value]= feval(Hfun, f,a,b,eps1,eps2,N);
8 E7 H) ^! v' f% p0 L
6 B, ]; [8 O$ ^/ w运行结果
2 y/ D6 ^0 V, v) O: Ek ,a ,b ,x, f
# C% H% Q5 W4 A4 o  1, 0.000000000, 1.000000000, 0.500000000,-352805.622314164,
' Y, {( k6 x4 i( W  ?  2, 0.000000000, 0.500000000, 0.250000000,-92968.331400711,
  3, 0.000000000, 0.250000000, 0.125000000,-4146.780462183,
0 O) f0 d$ ~8 O+ G  [  4, 0.000000000, 0.125000000, 0.062500000,29178.936016433,
5 o3 k+ `7 x" g9 `  5, 0.062500000, 0.125000000, 0.093750000,13458.364371693,
  6, 0.093750000, 0.125000000, 0.109375000,4890.124291139,
+ D" M7 }3 c! E' R* P  7, 0.109375000, 0.125000000, 0.117187500,430.092693351,
5 A% B# n6 @2 L% a( h  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,
$ _) W4 H7 z, R: z$ M4 y 11, 0.117187500, 0.118164063, 0.117675781,147.457778641,
, n) y2 |+ p# v9 z3 Y( v7 w# Z7 v. |* g2 h 12, 0.117675781, 0.118164063, 0.117919922,5.969190901,
13, 0.117919922, 0.118164063, 0.118041992,-64.817881139,
14, 0.117919922, 0.118041992, 0.117980957,-29.420780403,
15, 0.117919922, 0.117980957, 0.117950439,-11.724903562,
1 `9 Q' P2 q/ W& D 16, 0.117919922, 0.117950439, 0.117935181,-2.877633532,
+ A3 s2 f/ x; Z9 D8 p 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,
& z8 s: k8 z. \" @% _8 o7 ] 21, 0.117929459, 0.117930412, 0.117929935,0.163512627,
: P) C9 r. x# b# }2 A 22, 0.117929935, 0.117930412, 0.117930174,0.025279853,
23, 0.117930174, 0.117930412, 0.117930293,-0.043836575,
24, 0.117930174, 0.117930293, 0.117930233,-0.009278357,
25, 0.117930174, 0.117930233, 0.117930204,0.008000749,
: _0 T: K$ ~3 L, ~6 l2 l4 b: F5 f 26, 0.117930204, 0.117930233, 0.117930219,-0.000638804,
27, 0.117930204, 0.117930219, 0.117930211,0.003680972,
>>

8 ?2 A& }/ V1 F2 \- X; s结论 当x=-1.361时  y=0.117930219
* @( y0 W/ |9 J按此方法   可依次  得出150对x与y的值

单车居士

西交大学生？？问邢老师

明月山河

你这个是代数方程？

cosxuan

算法？

shouce

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

二分法        
* F7 M0 }7 z6 T% S还有一个程序
7 \/ {* t0 N! I" l$ k2 A4 \1 e+ d% Bisection.m
6 s* ~& G/ d9 {' m! kfunction[k,x,f_value]=Bisection(f,a,b,eps1,eps2,N)
fprintf('k ,a ,b ,x, f\n')
, K, m, A/ B6 efor 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)
: t( |3 F/ F3 X$ K  j5 ^    if abs(f_value)< eps1||0.5*(b-a)<eps2
$ Y2 J  v: [: S" P! a" \8 R    return
    else
  p' Q4 {% o1 f9 m        if f(x)*f(a)<0
; l3 J$ W, \( {2 C: h: A          b=x;
        else
+ y8 V. g, t5 D( N5 M$ d          a=x;
        end
        if k== N
          warning ('算法超出最大迭代数!')   
        end
0 r- s: P( c& M# F    end
; T0 }4 b! c% Z6 ]9 Y  D% D

shouce

单车居士 发表于 2015-6-4 12:41
4 V) c) Z; b! r1 U  }4 I  N$ e1 d西交大学生？？问邢老师

这点问题   自己就可以解决
4 Y8 g, z1 E2 j# [: c8 c# e$ L

shouce

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

是的呀      不过matlab写出来的方程  括号很多  不好看   

/ }2 M; t9 N9 l9 g

shouce

明月山河 发表于 2015-6-4 19:12
, ?. ]: b: I( C( b2 q# u2 ]你这个是代数方程？

' @3 P5 }5 B! B' x8 a
是超越方程      
& I2 i* C9 J7 ^3 V* F2 R1 k
  m0 B/ o3 ~0 V+ p
% c  H, a! N: f  ][chāo yuè fāng chéng
9 g$ L2 v4 I+ k; j  \: Y" Z. I: O- T