2011-10-22 4 views
0

Muller's method에 필요한 연산을 수행하는 루프가있는 함수를 작성하고 싶습니다.Mathematica의 뮬러 메서드

f[x_] := x^3 - x - 1; 
x0 = 0.8 
x1 = 1.5 
x2 = 2.0 
x3 = 5.0; 
\[Epsilon] = 0.001; 

While[(Abs[f[x3]] >= \[Epsilon]), 
h0 = x1 - x0; 
h1 = x2 - x1; 
d0 = (f[x1] - f[x0])/h0; 
d1 = (f[x2] - f[x1])/h1; 
A = (d1 - d0)/(h1 + h0); 
B = A*h1 + d1; 
Cx = f[x2]; 
raiz = Sqrt[B^2 - 4.0*A*Cx]; 
If[Abs[B + raiz] > Abs[B - raiz], dens = B + raiz, dens = B - raiz]; 
x3 = (x2 - 2*Cx)/dens; 
i++; 
Print["Iteration: ", i, "\t root \[TildeTilde] ", x3]; 
x0 = x1; 
x1 = x2; 
x2 = x3; 
] 

하지만 infinit 루프를 얻을 수는 ...

+0

'x3'은 (는)'x3 = x2-2 * Cx/dens'이어야한다고 생각합니다. – Heike

답변

4

뮬러의 방법 (항상 위키 백과보다 더) following Eric :

h[x_] := HermiteH[24, x]; 
i = [email protected][h[x], x] - 1; 
f[i, x_] := h[x]; 
roots = {}; 
While[ i > 1, 
    x0 = -2; x1 = -1; x2 = -.5; k = 1; 
    While[Abs[k] > .001, 
    q = (x0 - x1)/(x1 - x2); 
    a = q f[i, x0] - q (1 + q) f[i, x1] + q^2 f[i, x2]; 
    b = (2 q + 1) f[i, x0] - (1 + q)^2 f[i, x1] + q^2 f[i, x2]; 
    c = (1 + q) f[i, x0]; 
    p = Sqrt[b b - 4 a c]; 
    xp = x0 - (x0 - x1) 2 c /(k = If[Abs[b + p] > Abs[b - p], b + p, b - p]); 
    {x2, x1, x0} = {x1, x0, xp}; 
    ]; 
    AppendTo[roots, xp]; 
    i--; 
    f[i, x_] = f[i + 1, x]/(x - xp); 
    ]; 
Show[ 
Plot[h[x], {x, -2, 2}], 
Graphics[{PointSize[Large], Point[{#, 0} & /@ roots]}]] 
아래의 코멘트에 몇 가지 오류를 지적 헤이 케에 감사합니다

enter image description here

+0

이 방법은 수렴을 위해 세 가지 초기 점에 대한 좋은 추측을 요구합니다. –

+1

'k'의 정의가 맞다고 생각하지 않습니다. wikipedia 페이지와 Mathworld 페이지를 올바르게 해석하면'k'는'Abs [b + p]> Abs [b-p], b + p, b-p]'와 같아야합니다. – Heike

+0

나는 Abs [b + p]> Abs [b - p], k = b + p, k = b - p] 인 경우 제안을했으나 Divide :: infy : 무한 표현식 1/(0 . * 10^931 + 0. * 10^931 I). >>' – cMinor

관련 문제