Kolmiulotteinen Newtonin menetelmä

160 days ago by dx27

Newtonin iteraatio n-ulotteisessa tapauksessa:

x_{k+1}=x_k - J_f (x_k)^{-1}f(x_k)

Seuraavassa esimerkissä tarkastellaan funktiota

f(\mathbf{x})=\left(\begin{array}{c}3x_0-\cos{(x_1 x_2)}-\frac{1}{2}\\x_0^2-81(x_1+0.1)^2+\sin{x_2}+10.6\\e^{-x_0 x_1}+20x_2+\frac{10\pi-3}{3}\end{array}\right),

missä \mathbf{x}=(x_0, x_1, x_2). Ohjelmassa Jacobin matriisi lasketaan symbolisesti ja käänteismatriisi numeerisesti.

x0, x1, x2 = var('x0 x1 x2') f1(x0, x1, x2) = 3*x0 - cos(x1*x2) - (1/2) f2(x0, x1, x2) = x0^2 - 81*(x1 + 0.1)^2 + sin(x2) + 10.6 f3(x0, x1, x2) = e^(-x0*x1) + 20*x2 + (10*pi - 3)/3 f(x0, x1, x2) = (f1(x0,x1,x2), f2(x0,x1,x2), f3(x0,x1,x2)) j = jacobian(f, [x0,x1,x2]) x = vector([3.0, 4.0, 5.0]) data = [[0, x, n(norm(f(x[0], x[1], x[2])), digits=4)]] for i in range(1,8): x = vector((n(d) for d in x - j(x0=x[0], x1=x[1], x2=x[2]).inverse()*f(x[0], x[1], x[2]))) data.append([i, x, norm(f(x[0], x[1], x[2]))]) # Kuva viiva = line3d([d[1] for d in data], thickness=5) piste = point3d(data[-1][1], size=15, color="red") show(viiva + piste) # HTML-Taulukko s = "<br><br><table border=1 cellpadding=5>" s += "<tr><td> $ i $ </td><td> $ (x_0, x_1, x_2) $ </td><td> $ normi(f) $ </td></tr> " for i in range(0,8): s += "<tr><td> $ %s $ </td><td> $ %s $ </td><td> $ %s $ </td></tr>" %(data[i][0], data[i][1].n(digits=10), n(data[i][2], digits=4)) s += "</table>" html(s) 
       


i (x_0, x_1, x_2) normi(f)
0 (3.000000000, 4.000000000, 5.000000000) 1347.
1 (9.950035168, 2.038510442, -0.4735923501) 262.9
2 (0.5245370576, 0.7428920309, -0.4735987782) 47.13
3 (0.4949073626, 0.3972557032, -0.5143475312) 9.675
4 (0.4971709310, 0.2771459340, -0.5170792775) 1.169
5 (0.4970474904, 0.2580112360, -0.5175788308) 0.02966
6 (0.4970438097, 0.2574996297, -0.5175919449) 0.0001239
7 (0.4970438070, 0.2574992635, -0.5175919542) 0.0001221
\newcommand{\Bold}[1]{\mathbf{#1}}


i (x_0, x_1, x_2) normi(f)
0 (3.000000000, 4.000000000, 5.000000000) 1347.
1 (9.950035168, 2.038510442, -0.4735923501) 262.9
2 (0.5245370576, 0.7428920309, -0.4735987782) 47.13
3 (0.4949073626, 0.3972557032, -0.5143475312) 9.675
4 (0.4971709310, 0.2771459340, -0.5170792775) 1.169
5 (0.4970474904, 0.2580112360, -0.5175788308) 0.02966
6 (0.4970438097, 0.2574996297, -0.5175919449) 0.0001239
7 (0.4970438070, 0.2574992635, -0.5175919542) 0.0001221
\newcommand{\Bold}[1]{\mathbf{#1}}