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

shouce

螺杆转子刀具设计 数学模型

解如下方程：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
/ x6 b2 ~( i2 r# `: _1 Q求当x=-1.361时  y=?  
* D5 j7 Y" B6 j4 o) u! T>> syms x y
5 @( l* a3 r! E& [>> 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
0 s& p# m, J- h- x4 J
(x)+179999.4000005*y*sin(y)*sin(x)+75000*cos(y)*cos(x)*sin(x)

  N/ @2 M& S% `4 P" \, Ef =

3 p/ b) p2 B" y3 k(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)
  p. f+ t0 J" M
& B: A, t* c" c( m*sin(x) - (6184732290414159*y*cos(x)*cos(y))/34359738368 + (6184732290414159*y*sin(x)*sin(y))/34359738368

3 M3 p/ p- A* U1 V8 O>>  subs(f, x, -1.361)
; F) p; G8 ?3 d+ {5 S' L
ans =
  S* P  u: W' v) l
(3359995*sin(1361/1000)*cos(y))/2 - (2999995*cos(1361/1000)*sin(y))/2 - (7215543768789811*sin(1361/1000))/4294967296 + 75000*sin

7 z; L) Q5 l+ y, c6 r$ X, ](1361/1000)^3 - 75000*cos(1361/1000)*sin(1361/1000)*cos(y) - (6184732290414159*y*cos(1361/1000)*cos(y))/34359738368 -
6 Z2 u- B8 ?# f; b, q$ a8 [! \
3 a3 r  T/ @" I9 Q3 K: t0 Z3 o% |: x(6184732290414159*y*sin(1361/1000)*sin(y))/34359738368
* h, M# z1 j" q6 S' [% G1 u

) f) x% [' K9 c$ Y# c5 A>>  a=0; b=1;
eps1=1e-8;eps2=1e-8;
: L! T3 d, l- T% |) j" Q 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);
3 Q# m2 ]: I/ m: R  N8 A% ~Hfun=@Bisection;
[k,x,f_value]= feval(Hfun, f,a,b,eps1,eps2,N);
7 B: M, Z& d! c' K$ m( }6 E
- U2 G$ ^9 ]9 \, e5 w9 A9 G0 b- k运行结果
k ,a ,b ,x, f
  1, 0.000000000, 1.000000000, 0.500000000,-352805.622314164,
  2, 0.000000000, 0.500000000, 0.250000000,-92968.331400711,
: z3 J! G* W! V" E. x. z0 I  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,
( X9 O  `9 B& }  6, 0.093750000, 0.125000000, 0.109375000,4890.124291139,
  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,
; L. _# \6 C" N! O2 K 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,
13, 0.117919922, 0.118164063, 0.118041992,-64.817881139,
) u/ l1 O6 ~, F* d; s 14, 0.117919922, 0.118041992, 0.117980957,-29.420780403,
+ r: H; Z$ a$ |* S 15, 0.117919922, 0.117980957, 0.117950439,-11.724903562,
16, 0.117919922, 0.117950439, 0.117935181,-2.877633532,
: g4 h6 {" ^3 _$ j' Y# ` 17, 0.117919922, 0.117935181, 0.117927551,1.545834384,
18, 0.117927551, 0.117935181, 0.117931366,-0.665885649,
: H7 i! f3 c( T. Z$ k" U/ [5 q 19, 0.117927551, 0.117931366, 0.117929459,0.439977849,
/ |' ]; L( u) |  V4 e 20, 0.117929459, 0.117931366, 0.117930412,-0.112953030,
" G- M) D3 _. ~! Q! C2 f) W0 E& ^ 21, 0.117929459, 0.117930412, 0.117929935,0.163512627,
4 a" |* Z; S5 I/ ?6 ^7 f 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,
+ h+ d5 Z4 u6 m/ F 25, 0.117930174, 0.117930233, 0.117930204,0.008000749,
& H! n( Q" M1 k  [7 C/ K 26, 0.117930204, 0.117930233, 0.117930219,-0.000638804,
/ |* |! t; [, b7 b# M9 e 27, 0.117930204, 0.117930219, 0.117930211,0.003680972,
>>
0 {; b+ r$ h" `) z) G0 o: {
! B* G2 S- L* Q/ k结论 当x=-1.361时  y=0.117930219
按此方法   可依次  得出150对x与y的值
" Y( a) \9 d* u& f: p

单车居士

西交大学生？？问邢老师

明月山河

你这个是代数方程？

cosxuan

算法？

shouce

cosxuan 发表于 2015-6-4 19:45
1 P/ F( t" Y0 x1 [0 m' n8 P算法？

二分法        
还有一个程序
% Bisection.m
function[k,x,f_value]=Bisection(f,a,b,eps1,eps2,N)
- }1 W! j/ w' `- z( S( nfprintf('k ,a ,b ,x, f\n')
for k=1:N
    x=(a+b)/2;
    f_value=f(x);
    fprintf( '%3d, %10.9f, %10.9f, %10.9f,%10.9f,\n'...
1 O  ?  r5 y# x  q4 d      ,k ,a ,b ,x, f_value)
    if abs(f_value)< eps1||0.5*(b-a)<eps2
" t- F! z! R* P5 z1 h, _+ u    return
0 V* N/ k1 {$ T% c$ @( N% Z2 E    else
6 D# m3 t& F& y4 ?6 ~        if f(x)*f(a)<0
) J# _( b% K; u2 K$ o/ i9 f          b=x;
. T3 t7 k7 _' [& y) R        else
9 d) `$ g) O7 r* s$ d$ t          a=x;
        end
        if k== N
          warning ('算法超出最大迭代数!')   
        end
    end
/ {: S2 h; s2 l9 X" {; E% Z

shouce

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

这点问题   自己就可以解决
6 L2 F9 ?) p' {* K' C# A

shouce

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

% x! G" e- n' u, B3 R是的呀      不过matlab写出来的方程  括号很多  不好看   

shouce

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

% h& A+ K& ~  b6 m1 s) J$ _
是超越方程      


6 o, \; v: S2 h! ^8 ][chāo yuè fāng chéng
8 s+ u( `! m" {6 T