螺桿轉子刀具設計（2） matlab 解法

shouce

螺桿轉子刀具設計 數學模型
! U5 s& k6 ^1 Q4 D) y" L9 I
4 f4 ?! d; j$ Q# X, c解如下方程：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
$ T- C0 l. S) U$ v1 V求當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
8 M0 q) {% o- L, h+ M
(x)+179999.4000005*y*sin(y)*sin(x)+75000*cos(y)*cos(x)*sin(x)
3 g4 D# O7 Z8 I8 ~4 u
f =

(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)

: B5 n- [" ~1 d9 y8 R3 A; K*sin(x) - (6184732290414159*y*cos(x)*cos(y))/34359738368 + (6184732290414159*y*sin(x)*sin(y))/34359738368
4 S, j( x+ ]& E
>>  subs(f, x, -1.361)

% ~0 M; ~1 f7 K( w* a& u/ b# N2 xans =

' [9 ?( W& n8 D2 @8 n(3359995*sin(1361/1000)*cos(y))/2 - (2999995*cos(1361/1000)*sin(y))/2 - (7215543768789811*sin(1361/1000))/4294967296 + 75000*sin

. t# H; z: o; F" z  {: ](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
! f$ _3 ]/ Q* z3 |2 u3 J: r7 A
, X9 s8 u/ f$ q+ t! G
, S+ [9 c1 Y" w$ L>>  a=0; b=1;
! e7 d' e+ G% x' D1 Heps1=1e-8;eps2=1e-8;
+ {- W0 r( c' C1 b N=300;
9 P* n2 e8 ~3 @7 ?+ i+ i7 @9 f% Rf=@ (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);

運行結果
k ,a ,b ,x, f
  1, 0.000000000, 1.000000000, 0.500000000,-352805.622314164,
  2, 0.000000000, 0.500000000, 0.250000000,-92968.331400711,
" E1 N) A# N* m" }8 \" W0 s  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,
$ |1 E3 t* l  h6 b( U' x  7, 0.109375000, 0.125000000, 0.117187500,430.092693351,
7 L- A! J% P4 \! X) P8 Y" U  8, 0.117187500, 0.125000000, 0.121093750,-1843.759437944,
. S2 s. b+ H6 f& \/ V* w  9, 0.117187500, 0.121093750, 0.119140625,-703.184646645,
10, 0.117187500, 0.119140625, 0.118164063,-135.633470273,
6 y4 n2 K8 D0 D7 r& b1 @3 }- W 11, 0.117187500, 0.118164063, 0.117675781,147.457778641,
12, 0.117675781, 0.118164063, 0.117919922,5.969190901,
( O( x# Q% N2 I# {# o! P, V/ W 13, 0.117919922, 0.118164063, 0.118041992,-64.817881139,
9 c. W- ?* K; R; u& J9 c 14, 0.117919922, 0.118041992, 0.117980957,-29.420780403,
15, 0.117919922, 0.117980957, 0.117950439,-11.724903562,
9 O9 m# E! m# w, q0 s' J 16, 0.117919922, 0.117950439, 0.117935181,-2.877633532,
17, 0.117919922, 0.117935181, 0.117927551,1.545834384,
! o! m1 h) ~5 Q; H5 X 18, 0.117927551, 0.117935181, 0.117931366,-0.665885649,
  T) K- X1 Z4 B: b 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,
: X$ {- ?; l7 g( U" ^- T! w: ~ 23, 0.117930174, 0.117930412, 0.117930293,-0.043836575,
. R5 T7 Z4 S" o: |# h/ t 24, 0.117930174, 0.117930293, 0.117930233,-0.009278357,
' N% M* H/ Z# b: ^$ O  Z, D4 f) ~( M 25, 0.117930174, 0.117930233, 0.117930204,0.008000749,
26, 0.117930204, 0.117930233, 0.117930219,-0.000638804,
& j2 K) a( H& t$ Z  J: P 27, 0.117930204, 0.117930219, 0.117930211,0.003680972,
% ~$ s/ ?0 U& P9 r5 `' Z- ~>>

7 d; k) v; B- j, T' Z結論 當x=-1.361時  y=0.117930219
按此方法   可依次  得出150對x與y的值

單車居士

西交大學生？？問邢老師

明月山河

你這個是代數方程？

cosxuan

算法？

shouce

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

二分法        
9 X8 P$ v* T, |7 d- T3 ]還有一個程序
: h* u( k) _7 `" L+ B; {0 R5 `" Q! Q% Bisection.m
function[k,x,f_value]=Bisection(f,a,b,eps1,eps2,N)
fprintf('k ,a ,b ,x, f\n')
for k=1:N
5 t1 V  x( N2 U% e- \    x=(a+b)/2;
: D% Z! v" m+ F3 e    f_value=f(x);
% v3 [; V! K6 A, ?4 b9 ^    fprintf( '%3d, %10.9f, %10.9f, %10.9f,%10.9f,\n'...
2 K, c* {1 E, ~- J) ]      ,k ,a ,b ,x, f_value)
2 y6 G& A: R' ?1 y    if abs(f_value)< eps1||0.5*(b-a)<eps2
    return
; v# p* e  }+ P  F    else
" d8 W/ k  d$ @' Z$ h        if f(x)*f(a)<0
: m. J8 s& l  x  e- Y          b=x;
( W+ k$ `. [7 [/ b- ^3 ]. d        else
          a=x;
* ^: a! F; b- _! Q& a+ u2 i        end
        if k== N
          warning ('算法超出最大迭代數!')   
3 \# |: p1 c, W$ B! d        end
    end

shouce

單車居士 發表于 2015-6-4 12:41
0 i2 F3 S: r2 s- S: \: T+ J9 ?) ]西交大學生？？問邢老師

  J4 L9 e6 U' j這點問題   自己就可以解決

shouce

明月山河 發表于 2015-6-4 19:12
你這個是代數方程？

* H1 t0 ?1 b4 O: C& x0 a  m. X$ u3 ^是的呀      不過matlab寫出來的方程  括號很多  不好看   
9 A5 V- a' K) t  @
+ ]9 \8 `  y1 X5 b9 J' r

shouce

明月山河 發表于 2015-6-4 19:12
0 |% o1 K# d' Z, r/ v你這個是代數方程？

3 s5 i% B; {6 a& c" p$ r
7 M* i% ~' s9 a* N% G是超越方程      
( x) S! a: S* y
; D% e2 C* @! `& ?+ A5 t
[chāo yuè fāng chéng
& X$ f5 `, m2 Z4 o