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

shouce

螺杆转子刀具设计 数学模型
+ x1 E! N$ O' e, a5 w* |1 X/ 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
>> 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

4 ^- E% q. r) \3 Y! V$ o(x)+179999.4000005*y*sin(y)*sin(x)+75000*cos(y)*cos(x)*sin(x)
8 m# m5 F9 T! _* H
f =

/ u% O# F$ Y4 r1 M$ Y(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)

*sin(x) - (6184732290414159*y*cos(x)*cos(y))/34359738368 + (6184732290414159*y*sin(x)*sin(y))/34359738368

* B+ m8 \5 O- R( R7 n5 X2 ?>>  subs(f, x, -1.361)
; x4 @* o, s6 W3 M0 n
9 r( I7 C) Z7 B% {ans =
' f$ ^1 i# E, ^( }9 _' c
* J5 E+ w" y$ }# ]& f& g" b  E2 W(3359995*sin(1361/1000)*cos(y))/2 - (2999995*cos(1361/1000)*sin(y))/2 - (7215543768789811*sin(1361/1000))/4294967296 + 75000*sin
" F) L" e) l8 t% R
(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
& b. T* q( r( A+ @, }, f
' a4 z# t6 v9 v: L; I3 n( s
" ?! C' M0 A2 V9 c# z. j>>  a=0; b=1;
5 r, H$ ^: r' O0 n- u& ^% o, Ceps1=1e-8;eps2=1e-8;
$ `) K3 Q  h6 K 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);
% p/ f5 x+ C1 E% nHfun=@Bisection;
[k,x,f_value]= feval(Hfun, f,a,b,eps1,eps2,N);
+ c1 o, p- H, \3 m& A* V
运行结果
k ,a ,b ,x, f
  1, 0.000000000, 1.000000000, 0.500000000,-352805.622314164,
( \- Z% ^& u* m) N; j; f% t  2, 0.000000000, 0.500000000, 0.250000000,-92968.331400711,
  3, 0.000000000, 0.250000000, 0.125000000,-4146.780462183,
% o/ C; J: v4 P! {  F0 g  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,
1 `: p0 p; [& z/ }* _1 g  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,
, @" |8 l2 l4 o- Z 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,
14, 0.117919922, 0.118041992, 0.117980957,-29.420780403,
( g$ V8 o; G$ X( c4 j. Y& U 15, 0.117919922, 0.117980957, 0.117950439,-11.724903562,
, H- G, Y7 j, J 16, 0.117919922, 0.117950439, 0.117935181,-2.877633532,
- F  T+ c% O4 f/ J7 x+ V; s8 Q 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,
21, 0.117929459, 0.117930412, 0.117929935,0.163512627,
22, 0.117929935, 0.117930412, 0.117930174,0.025279853,
; ~; n% H2 Y* L& J0 q  x( x 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,
; o4 S7 N( |7 ^% i/ B+ J( X 26, 0.117930204, 0.117930233, 0.117930219,-0.000638804,
27, 0.117930204, 0.117930219, 0.117930211,0.003680972,
>>

9 M# G3 H& d( u: l2 n结论 当x=-1.361时  y=0.117930219
按此方法   可依次  得出150对x与y的值
% ~( C7 N# b0 p

单车居士

西交大学生？？问邢老师

明月山河

你这个是代数方程？

cosxuan

算法？

shouce

cosxuan 发表于 2015-6-4 19:45
3 ?$ {) Y; [+ f* v, M# t算法？

二分法        
6 k% x- D5 P6 ?$ X% o还有一个程序
% Bisection.m
' z7 \7 n' A. L# `/ lfunction[k,x,f_value]=Bisection(f,a,b,eps1,eps2,N)
# G5 `. s$ y5 r0 L) O& L. k/ a2 |fprintf('k ,a ,b ,x, f\n')
for k=1:N
    x=(a+b)/2;
( U, P1 m' y( S# V0 }    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
# i5 w" i' a# Y9 B0 T# e& D    return
    else
        if f(x)*f(a)<0
          b=x;
9 R  l- R$ j! m; T" Z! U/ Q        else
3 l/ p7 u! P1 p' ~          a=x;
! U/ O1 c) o1 x0 ~0 L        end
' W% B9 B  j" V" {7 R        if k== N
          warning ('算法超出最大迭代数!')   
        end
7 B6 R& e: H  _+ \    end

shouce

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

8 X  `) z9 g9 }- M/ k6 ~; [这点问题   自己就可以解决
' A1 H2 F1 E) E

shouce

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

是的呀      不过matlab写出来的方程  括号很多  不好看   
9 M$ X8 o- g% S) P# ?# c5 S' ^

shouce

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

' {9 `0 L+ N: }# o1 J! x
6 \7 ?0 b0 i1 ^/ N, R是超越方程      

2 H$ |5 e2 Q+ a4 b& v
[chāo yuè fāng chéng