function [x, f, ok, count,x2] = newton1D(fun, funprim, funbis, xinit, epsi, maxiter)

x2=[];

if nargin < 6, maxiter = 30; end;

if nargin < 5, epsi = 100 * eps; end;

x0 = xinit;

count = 0;

while 1

if count > maxiter

break

end;

fbis = feval(funbis, x0);

if fbis == 0

ok = 0;

disp('zerowa druga pochodna w newton1D');

return

end;

fprim = feval(funprim, x0);

x = x0 - fprim / fbis;

x2=[x2 x];

if abs(x - x0) <= epsi * abs(x0)

if fbis > 0

ok = 1;

else

ok = 2;

end;

f = feval(fun, x);

return

end;

x0 = x;

count = count + 1;

end;

warning('przekroczona liczba iteracji w newton1D');

ok = 0;

Zadanie 1

Z zastosowaniem procedury newton1D określić minimum funkcji f(x)=x²/2 - sinx z dokładnością do czterech miejsc po przecinku. Jako punkt startowy wybrać x=0.5. co otrzymuje się po zmianie punku startowego? Narysować wykresy charakteryzujące zbieżność metody ( wartości x i f(x) względem liczby wykonanych iteracji) wykonać to samo dla funkcji y=sinx.

W podobny sposób przeanalizować pracę algorytmu dla funkcji y=x arctg(x)-0.5ln(1+x²),

badając szczególnie wrażliwość na zmianę punktu startowego. Wyjaśnić dziwne zachowanie metody gdy np.x^0=1.5.

Dla funkcji f(x)=x²/2 - sinx

script_zad1.

fun = inline('(x.^2)/2-sin(x)');

funprim = inline('x-cos(x)');

funbis = inline('1+sin(x)');

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 0.5, 1.e-6, 100)

script_zad1;

x =

0.7391

f =

-0.4005

ok =

1

count =

3

x2 =

0.7552 0.7391 0.7391 0.7391

y2=fun(x2);

y2 =

-0.4003 -0.4005 -0.4005 -0.4005

plot(x2,y2,'+');

hold on

plot(x2,y2,'+');

x3=0.5:0.001:1;

y3=fun(x3);

plot(x3,y3);

Wykres 1

script_zad1.

fun = inline('(x.^2)/2-sin(x)');

funprim = inline('x-cos(x)');

funbis = inline('1+sin(x)');

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 1.5, 1.e-6, 100)

script_zad1_1;

x =

0.7391

f =

-0.4005

ok =

1

count =

3

x2 =

0.7845 0.7395 0.7391 0.7391

y2=fun(x2)

y2 =

-0.3988 -0.4005 -0.4005 -0.4005

plot(x2,y2,'+');

hold on

x2=0.5:0.01:1;

y2=fun(x2);

plot(x2,y2);

Dla funkcji y=sinx.

script_zad1_2

fun = inline('sin(x)');

funprim = inline('cos(x)');

funbis = inline('-sin(x)');

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 0.5, 1.e-6, 100)

script_zad1_2;

x =

1.5708

f =

1

ok =

2

count =

4

x2 =

2.3305 1.3806 1.5731 1.5708 1.5708

y2=feval(fun,x2)

y2 =

0.7250 0.9820 1.0000 1.0000 1.0000

plot(x2,y2,'+')

hold on

plot(x3,y3);

x3=1:0.01:3

y3=fun(x3);

plot(x3,y3);

Wykres 2

script_zad1_2;

x =

1.5708

f =

1

ok =

2

count =

2

x2 =

1.5709 1.5708 1.5708

y2=fun(x2)

y2 =

1.0000 1.0000 1.0000

plot(x2,y2,'+');

hold on

x3=1.5:0.001:1.6

y3=fun(x3);

plot(x2,y2,'+');

hold on

plot(x3,y3);

Dla funkcji y=x arctg(x)-0.5ln(1+x˛)

script_zad1_3

fun = inline('x.*atan(x)-0.5*log(1+x.^2)');

funprim = inline('atan(x)');

funbis = inline('1/(1+x.^2)');

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 0.5, 1.e-6, 100)

script_zad1_3

x =

0

f =

0

ok =

1

count =

4

x2 =

-0.0796 0.0003 0.0000 0 0

y2=fun(x2);

plot(x2,y2,'+');

hold on;

x3=-0.09:0.001:0.01;

y3=fun(x3);

plot(x3,y3);

Wykres 3

Dla tej samej funkcji ale z innym punktem startowym. X⁰=1.5

x =

-9.4595e+216

ok =

0

count =

11

x2 =

1.0e+216 *

Columns 1 through 7

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

Columns 8 through 11

0.0000 0.0000 0.0000 -9.4595

y2 = fun(x2);

y2 =

1.0e+108 *

plot(x2,y2,'+');

hold on;

x3=-1.5:0.001:2;

y3=fun(x3);

plot(x3,y3);

x3=-1.5:0.001:2

Wykres

Zadanie 2 (dla fanów)

Funkcja

f(x)= - + - +

posiada cztery ekstrema w przedziale [0;1]; zauważyć, ze

f'(x) = (x-1/2) * (x-1/3) * (x-1/4) * (x-1/5)

Przedstawić na wykresie zależność punktu ekstremum, do którego jest zbieżna metoda Newtona, od punktu startowego ( założyć, że punkty startowe wybiera się z przedziału [0;1].

Dla funkcji f(x)= - + - +

fun = inline('(x.^5)/5-(77.*x.^4)/240+(71.*x.^3)/360-(7.*x.^2)/120+(x/120)');

funprim = inline('(x-1/2).*(x-1/3).*(x-1/4).*(x-1/5)');

funbis = inline('(x-1/3).*(x-1/4).*(x-1/5)+(x-1/2).*(x-1/4).*(x-1/5)+(x-1/2).*(x-1/3).*(x-1/5)+(x-1/2).*(x-1.3).*(x-1/4)');

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 0.0, 1.e-6, 100)

y2=fun(x2)

plot(x2,y2,'+m');

hold on;

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 0.2, 1.e-6, 100)

y2=fun(x2)

plot(x2,y2,'+r');

hold on;

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 0.4, 1.e-6, 100)

y2=fun(x2)

plot(x2,y2,'+g');

hold on;

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 0.6, 1.e-6, 100)

y2=fun(x2)

plot(x2,y2,'+y');

hold on;

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 0.8, 1.e-6, 100)

y2=fun(x2)

plot(x2,y2,'+b');

hold on;

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 1, 1.e-6, 100)

y2=fun(x2)

plot(x2,y2,'+c');

hold on;

x3=0:0.01:1;

y3=fun(x3);

plot(x3,y3);

script_zad2

x =

0.2000

f =

4.6178e-004

ok =

2

count =

87

x2 =

Columns 1 through 7

0.0351 0.0630 0.0853 0.1034 0.1182 0.1304 0.1405

Columns 8 through 14

0.1490 0.1561 0.1621 0.1672 0.1716 0.1753 0.1785

Columns 15 through 21

0.1813 0.1837 0.1858 0.1876 0.1891 0.1905 0.1917

Columns 22 through 28

0.1927 0.1936 0.1944 0.1951 0.1957 0.1962 0.1967

Columns 29 through 35

0.1971 0.1974 0.1978 0.1980 0.1983 0.1985 0.1987

Columns 36 through 42

0.1988 0.1990 0.1991 0.1992 0.1993 0.1994 0.1995

Columns 43 through 49

0.1995 0.1996 0.1996 0.1997 0.1997 0.1998 0.1998

Columns 50 through 56

0.1998 0.1998 0.1999 0.1999 0.1999 0.1999 0.1999

Columns 57 through 63

0.1999 0.1999 0.1999 0.1999 0.2000 0.2000 0.2000

Columns 64 through 70

0.2000 0.2000 0.2000 0.2000 0.2000 0.2000 0.2000

Columns 71 through 77

0.2000 0.2000 0.2000 0.2000 0.2000 0.2000 0.2000

Columns 78 through 84

0.2000 0.2000 0.2000 0.2000 0.2000 0.2000 0.2000

Columns 85 through 88

0.2000 0.2000 0.2000 0.2000

y2 =

1.0e-003 *

Columns 1 through 7

0.2286 0.3378 0.3928 0.4218 0.4377 0.4469 0.4523

Columns 8 through 14

0.4556 0.4576 0.4590 0.4598 0.4604 0.4608 0.4611

Columns 15 through 21

0.4613 0.4614 0.4615 0.4616 0.4616 0.4617 0.4617

Columns 22 through 28

0.4617 0.4617 0.4617 0.4618 0.4618 0.4618 0.4618

Columns 29 through 35

0.4618 0.4618 0.4618 0.4618 0.4618 0.4618 0.4618

Columns 36 through 42

0.4618 0.4618 0.4618 0.4618 0.4618 0.4618 0.4618

Columns 43 through 49

0.4618 0.4618 0.4618 0.4618 0.4618 0.4618 0.4618

Columns 50 through 56

0.4618 0.4618 0.4618 0.4618 0.4618 0.4618 0.4618

Columns 57 through 63

0.4618 0.4618 0.4618 0.4618 0.4618 0.4618 0.4618

Columns 64 through 70

0.4618 0.4618 0.4618 0.4618 0.4618 0.4618 0.4618

Columns 71 through 77

0.4618 0.4618 0.4618 0.4618 0.4618 0.4618 0.4618

Columns 78 through 84

0.4618 0.4618 0.4618 0.4618 0.4618 0.4618 0.4618

Columns 85 through 88

0.4618 0.4618 0.4618 0.4618

x =

0.2000

f =

4.6178e-004

ok =

2

count =

0

x2 =

0.2000

y2 =

4.6178e-004

x =

0.5000

f =

4.3403e-004

ok =

1

count =

7

x2 =

Columns 1 through 7

0.4179 0.4416 0.4698 0.4926 0.4996 0.5000 0.5000

Column 8

0.5000

y2 =

1.0e-003 *

Columns 1 through 7

0.4528 0.4462 0.4384 0.4343 0.4340 0.4340 0.4340

Column 8

0.4340

x =

0.5000

f =

4.3403e-004

ok =

1

count =

3

x2 =

0.5004 0.5000 0.5000 0.5000

y2 =

1.0e-003 *

0.4340 0.4340 0.4340 0.4340

x =

0.5000

f =

4.3403e-004

ok =

1

count =

4

x2 =

0.6185 0.5067 0.4998 0.5000 0.5000

y2 =

1.0e-003 *

0.6541 0.4343 0.4340 0.4340 0.4340

x =

0.5000

f =

4.3403e-004

ok =

1

count =

5

x2 =

0.7659 0.5943 0.4989 0.5000 0.5000 0.5000

y2 =

0.0031 0.0006 0.0004 0.0004 0.0004 0.0004

Wykres

Zadnie 3

Określić punkty minimum i maksimum funkcji

1. f(x)=x(x-1)²

2. f(x)=x/(x ²+1)

a) Maksimum dla funkcji f(x)=x(x-1)²

fun = inline('x.*(x-1).^2');

funprim = inline('(x-1).^2+2.*x.*(x-1)');

funbis = inline('4.*(x-1)+2.*x');

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 0, 1.e-6, 100)

zad3

x =

0.3333

f =

0.1481

ok =

2

count =

4

x2 =

0.2500 0.3250 0.3332 0.3333 0.3333

y2=fun(x2);

plot(x2,y2,'+');

hold on

x3=0.2:0.001:0.5;

y3=fun(x3);

plot(x3,y3)

Wykres

Minimum dla funkcji f(x)=x/(x ˛+1)

fun = inline('x.*(x-1).^2');

funprim = inline('(x-1).^2+2.*x.*(x-1)');

funbis = inline('4*(x-1)+2*x');

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 1.5, 1.e-6, 100)

d 1.5

|

Missing operator, comma, or semi-colon.

zad3

x =

1.0000

f =

1.2509e-024

ok =

1

count =

4

x2 =

1.1500 1.0233 1.0008 1.0000 1.0000

y2=fun(x2)

y2 =

0.0259 0.0006 0.0000 0.0000 0.0000

plot(x2,y2,'+');

plot(x2,y2,'+');

x3=0.7:0.01:1.3

y3=fun(x3)

hold on

plot(x3,y3)

Wykres

b)

fun = inline('x./(x.^2+1)');

funprim = inline('(-x.^2+1)/(x.^2+1).^2');

funbis = inline('(2.*x.^5-4.*x.^3-6.*x)/(x.^2+1).^4');

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 0.5, 1.e-6, 100)

zad3b

x =

1

f =

0.5000

ok =

2

count =

5

x2 =

0.8409 0.9706 0.9988 1.0000 1.0000 1.0000

y2=fun(x2)

y2 =

0.4926 0.4998 0.5000 0.5000 0.5000 0.5000

plot(x2,y2,'+');

x3=0.8:0.01:1.2;

y3=fun(x3);

hold on

plot(x3,y3);

fun = inline('x./(x.^2+1)');

funprim = inline('(-x.^2+1)/(x.^2+1).^2');

funbis = inline('(2.*x.^5-4.*x.^3-6.*x)/(x.^2+1).^4');

[x, f, ok, count,x2] = newton1D(fun, funprim, funbis, 1.5, 1.e-6, 100)

Dla 1.5

dla 1.5

|

Missing operator, comma, or semi-colon.

zad3b

x =

-1.0000

f =

-0.5000

ok =

1

count =

5

x2 =

-0.3056 -0.8636 -0.9776 -0.9993 -1.0000 -1.0000

y2=fun(x2)

y2 =

-0.2795 -0.4947 -0.4999 -0.5000 -0.5000 -0.5000

x3=-1.2:0.01:0;

y3=fun(x3);

plot(x2,y2,'+');

hold on

plot(x3,y3)

---