Maplet Programming

Phương pháp Newton cho hệ phi tuyến


Lý thuyết xấp xỉ nghiệm của hệ phi tuyến tổng quát bằng phương pháp Newton có thể xem trong bài báo dưới đây 

http://benisrael.net/Newton-MP.pdf

Ở đây giới thiệu maplet cho bài xấp xỉ nghiệm của hệ phi tuyến bẳng phương pháp Newton hiểu theo nghĩa cổ điển.

Thí dụ giải hệ phương trình x^3-y+1=0, 2x+y+3=0

Trước hết lập một vài procedure cần thiết

> with(LinearAlgebra): with(MTM):
> F:=proc(F1,F2,x,y) local f1,f2;      # Lập vector hàm
f1:= proc (x,y) options operator, arrow; F1 end proc;
f2:= proc (x,y) options operator, arrow; F2 end proc;
return (Vector([f1(x,y),f2(x,y)]));
end proc:

> JacobianMatrix := proc (F1,F2,x,y) local f1,f2;    # tính ma trận Jacobian
f1:= proc (x,y) options operator, arrow; F1 end proc;
f2:= proc (x,y) options operator, arrow; F2 end proc;
return (jacobian(Vector([f1(x,y),f2(x,y)]),Vector([x,y])))
end proc:

> V := proc (n::nonnegint,F1,F2,x,y)  #  Tìm nghiệm gần đúng bằng đệ quy
if n = 1 then return Vector([-1, 1]);             
         else return (
V(n-1,F1,F2,x,y)-MatrixVectorMultiply(
MatrixInverse(eval(JacobianMatrix(F1,F2,x,y),[x=V(n-1,F1,F2,x,y)[1],y=V(n-1,F1,F2,x,y)[2]])),
eval(F(F1,F2, x,y),[x=V(n-1,F1,F2,x,y)[1],y=V(n-1,F1,F2,x,y)[2]] ))                     );
end if;
end proc:

> U:=proc(F1,F2,x,y,e) local n;  # Tìm nghiệm gần đúng với sai số là  e
n:=1;
while (VectorNorm(V(n+1,F1,F2,x,y)-V(n,F1,F2,x,y))> e) do
           n:=n+1; end do;  
return(evalf(V(n,F1,F2,x,y)));
end proc:

Có thể test một số hệ phương trình hai biến (dĩ nhiên có thể xảy ra trường hợp phương trình vô nghiệm, hoặc không hội tụ, tuy nhiên trong giải thuật ở trên được xin bỏ qua những điều này, việc bổ sung và tổng quát giải thuật không có nhiều khó khăn), chẳng hạn với ví trên

>  U(x^3-y+1, 2*x+y+3, x, y, 1/1000000);

Kết quả:  x=-1.179509057, y=-0.6409818860

Lập một maplet là bước cuối cùng để thu được sản phẩm phần mềm hữu ích
> with(Maplets[Elements]):

> maplet := Maplet(
Window(‘title’ = “Solving System of Equations”, ‘toolbar’ = ‘TB’,
[[“Enter the 1st function with variables x,y”], [TextField[‘f1’](30)],
[“Enter the 2nd function with variables x,y”], [TextField[‘f2’](30)],
[“Enter tolerance”, TextField[‘e’](12)],
[“Solution”],[MathMLViewer[‘ML’]()]]),

ToolBar[TB](
ToolBarButton(caption = “Do It”, image = Image(“run.png”), onclick=Evaluate(‘ML’=’U(f1,f2,x,y,e)’)),
ToolBarSeparator(),
ToolBarButton(caption = “About”, image = Image(“about.png”), onclick = RunDialog(‘DG’)),
ToolBarSeparator(),
ToolBarButton(caption = “Exit”, image = Image(“exit.png”), onclick = Shutdown())),
 
MessageDialog[DG](“Solving System of Equations. Copyringt (C) 2009 Mathematical Diary”)):

> Maplets[Display](maplet);

systemeq