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

2009 October 7
1 Comment
tags: , ,
by Tuan Minh

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

from → Maplet Programming

Circulant Matrices

2009 July 17
Leave a comment
tags: ,
by Tuan Minh

Definition.  Given T:\mathbb{C}^n\to \mathbb{C}^n be a linear operator such that for every x=(x_0,x_2,...,x_{n-1}) then

T(x_0,x_1,...,x_{n-1})=(x_{n-1},x_0,x_1,...,x_{n-2})

The matrix with k-row is given by T^{k-1}x, \ k=1,2,...,n will be called circulant and denoted by X=circ(x)=circ(x_0,x_1,...,x_n)

X=\begin{pmatrix}x_0 & x_1 & ... & x_{n-2} & x_{n-1}\\ x_{n-1} & x_0 & ... & x_{n-3} & x_{n-2}\\ ... & ... & ... & ... & ...\\ x_2 & x_3 &... & x_0 & x_1\\ x_1 & x_2 & ... & x_{n-1} & x_0\end{pmatrix}

Theorem 1. \det(X)=\prod_{k=0}^{n-1}\left ( \sum_{j=0}^{n-1}x_j \epsilon_k^j\right ), where \epsilon is a primitive n-th root of unity.

We multiply the  j-columm with \epsilon^{j-1}, \ j=2,3,...,n, and add it onto the first columm, claim that \epsilon^n=1, we have:

 

The j-element of the first columm is a product of x_0+x_1\epsilon+...+x_{n-1}\epsilon^{n-1} with \epsilon^{j-1}, \ j=1,2,...,n. We give \det(X) as a polynomial with n variable x_1,x_2,...,x_n on \mathbb{C}[x], hence \det(X) is devided  by x_0+x_1\epsilon+...+x_{n-1}\epsilon^{n-1}. Claim that \det(X) is a n-degree polynomial and the coefficient of x_0x_1...x_{n-1} is 1.

Thus \det(X)=\prod_{k=0}^{n-1}\left ( \sum_{j=0}^{n-1} \epsilon_k^jx_j \right ) \bigstar

For k=0,1,2,...,n-1, we denote  e_k=\frac{1}{\sqrt{n}}(1,\epsilon_k,...,\epsilon_k^{n-1})^T and
\lambda_k= x_0+x_1\epsilon_k+...+x_{n-1}\epsilon_k^{n-1}

We can easily calculate and receive that Xe_k=\lambda_k e_k and since \{e_0,e_1,...,e_{n-1}\} is a linearly independent system, we conclude that \lambda_0,\lambda_1,...,\lambda_{n-1} are eigenvalues of circulant matrix X

Some results:

i.  trace(X)=\sum_{k=0}^{n-1}\lambda_k=nx_0

ii. With non-singular matrix

 P=\frac{1}{\sqrt{n}} \begin{pmatrix}1 & 1 & ... & 1 & 1\\ 1 & \epsilon & ... & \epsilon^{n-2} & \epsilon^{n-1}\\ ... & ... & ... & ... & ...\\ 1 & \epsilon^{n-2} & ....& \epsilon^{(n-2)^2}& \epsilon^{(n-1)(n-2)}\\ 1 & \epsilon^{n-1} & ... & \epsilon^{(n-2)(n-1)} & \epsilon^{(n-1)^2} \end{pmatrix}

then P^{-1}XP=diag(\lambda_0,\lambda_1,...,\lambda_{n-1})

iii  X=\sum_{k=0}^{n-1}x_kW^k where W=circ(0,1,0,...,0)

(To be continued)

from → Matrix Theory

Chéo hóa ma trận bằng phương pháp Jacobi

2009 July 15
Leave a comment
tags:
by Tuan Minh

Có thể nói đây là phương pháp đơn giản nhất đối với bài toán chéo hóa các ma trận đối xứng.

Giả sử ma trận A\in M_n(\mathbb{R}) đối xứng. Ta xây dựng dãy ma trận trực giao (R_k)_{k\ge1} và dãy ma trận (D_k)_{k\ge 0} như sau

 

D_0=A, \ D_k=R_k^T D_{k-1}R_k, \ \ k=1,2,...

 

sao cho \displaystyle\lim_{k\to\infty} D_k = diag(\lambda_1,\lambda_2,...,\lambda_n) là dạng chéo hóa của A.

 

Thiết lập một dãy như vậy được gọi là quá trình lặp Jacobi

 

Tại bước thứ k của quá trình này D_{k}=R_{k}^T D_{k-1}R_k chúng ta sẽ làm cho phần tử d^{(k)}_{pq}=d^{(k)}_{qp} nằm ngoài đường chéo chính nhận giá trị 0 bằng cách đặt

 

R_k=\begin{pmatrix}1 & 0 & ... & 0 & ... & 0 &... & 0 & 0\\ 0 & 1 & ... & 0 & ... & 0 &... & 0 & 0\\ ...& ... & ... & ... & ... & ... &... & ... & ...\\ 0 & 0 & ... & c & ... & s &... & 0 & 0\\ ...& ... & ... & ... & ... & ... &... & ... & ...\\ 0 & 0 & ... & -s & ... & c &... & 0 & 0\\ ...& ... & ... & ... & ... & ... &... & ... & ...\\ 0 & 0 & ... & 0 & ... & 0 &... & 1 & 0\\ 0 & 0 & ... & 0 & ... & 0 &... & 0 & 1\end{pmatrix}

 

trong đó c=\cos\varphi, s=\sin\varphi

 

Tương ứng với phép biến đổi

d^{(k)}_{jp}=cd^{(k-1)}_{jp}-sd^{(k-1)}_{jq}, \ j\neq p, \ j\neq q

d^{(k)}_{jq}=s d^{(k-1)}+c d^{(k-1)}_{jq},  j\neq p, \ j\neq q

d^{(k)}_{pp}=c^2d^{(k-1)}_{pp}+s^2 d^{(k-1)}_{qq}-2csd^{(k-1)}_{pq}

 

d^{(k)}_{qq}=c^2d^{(k-1)}_{pp}+s^2 d^{(k-1)}_{qq}+2csd^{(k-1)}_{pq}

 

d^{(k)}_{pq}=(c^2-s^2)d^{(k-1)}_{pq}+cs\left( d^{(k-1)}_{pp}-d^{(k-1)}_{qq} \right)

 

Giả sử d^{(k-1)}_{pq}\neq 0, chúng ta muốn rằng d^{(k)}_{pq}=0 nên từ phương trình cuối cùng ta có

 

\theta=\cot2\varphi=\frac{c^2-s^2}{cs}=\frac{ d^{(k-1)}_{pp}-d^{(k-1)}_{qq}}{d^{(k-1)}_{pq}}

 

Ta đặt \displaystyle t=\tan\varphi=\frac{s}{c} suy ra \displaystyle\theta=\frac{1-s^2/c^2}{2s/c}=\frac{1-t^2}{2t}, dẫn đến phương trình bậc 2

 

t^2+2\theta t-1=0

 

với nghiệm \displaystyle t=\frac{sign(\theta)}{|\theta|+\sqrt{\theta^2+1}}. Từ đây ta tính được

 

c=\displaystyle\frac{1}{\sqrt{1+t^2}}, \ \ s=ct

 

Chúng ta xét các tổng các phần tử nằm ngoài đường chéo chính

 

S_k= \sum_{i\neq j, \ i,j=1...n} |d^{(k)}_{ij}|^2

 

Khi đó theo hệ phương trình đã xét ta có S_k=S_{k-1}-2|d^{(k-1)}_{ij}|^2

 

Như vậy dãy (S_k)_{k\ge 1} đơn điệu giảm và bị chặn dưới bởi 0 nên nó hội tụ.

 

Với vị trí (p,q) cần đưa về d^{(k)}_{pq}=0 trong ma trận D_k ta chọn tại vị trí của phần tử \displaystyle\max_{p<q}\{ |d^{(k-1)}_{pq}|\} trong ma trận D_{k-1}.

 

Với cách chọn như vậy thì (S_k)_{k\ge 1} hội tụ về 0 và hệ quả là (D_k)_{k\ge 0} hội tụ về ma trận đường chéo K. Ta biết rằng đa thức đặc trưng ma trận luôn bất biến với phép đồng dạng nên K chính là dạng chéo hóa của A.

 

Quá trình này được cài đặt với Mathematica như sau

 

Jacobi1 

 Jacobi2

Ví dụ.  

A=\begin{pmatrix} 8 & -1 & 3 & -1 \\ -1 & 6 & 2 & 0 \\ 3 & 2 & 9 & 1 \\ -1 & 0 & 1 & 7 \end{pmatrix}

 

Ta làm như sau:

 

{K, V} = JacobiCyclic[A, 10^-11, 50];

Print[" K=", MatrixForm[Chop[K]]];

Print[" V=", MatrixForm[Chop[V]]];

 

Khi đó kết quả thu được là dạng chéo hóa

 

K=\begin{pmatrix} 6.59234 & 0 & 0 & 0 \\ 0 & 3.2957 & 0 & 0 \\ 0 & 0 & 11.7043 & 0 \\ 0 & 0 & 0 & 8.40766 \end{pmatrix}

 

và ma trận trực giao

 

V=\begin{pmatrix} 0.230097 & 0.528779 & 0.582298 & -0.573042 \\ 0.628975 & 0.591967 & 0.175776 & 0.472301 \\ -0.0712347 & -0.536039 & 0.792487 & 0.28205 \\ 0.739169 & 0.287455 & 0.0446803 & 0.607455\end{pmatrix}

 

thỏa mãn A=V^T K V.

 

from → Matrix Theory

Popup Menu

2009 July 5
Leave a comment
tags: , ,
by Tuan Minh

PopupMenu là một đối tượng dạng menu xuất hiện khi ta kích chuột phải trên cửa sổ ứng dụng, hoặc cục bộ trên các đối tượng TextBox hoặc TextField. Như vậy nó là được coi như là một option của Window hoặc của đối tượng TextBox, TextField. PopupMenu bao gồm các Menu, MenuSeparator, MenuItem.

Thí dụ sau đây về ứng dụng tính đạo hàm của hàm hai biến. Trong đó một TextField để input hàm hai biến cần tính, PopupMenu của nó gồm:

- Menu: gồm hai ItemMenu tương ứng với hai chứng năng là tính đạo hàm theo x hoặc theo y và output ra một MathMLViewer.

- MenuSeparator (dấu ngăn cách)

- MenuItem: Thoát khỏi ứng dụng.

Thực hiện như sau

> with(Maplets[Elements]):

> MyMaplet := Maplet(Window(title=”Differentiate Calculator”,

[

["Enter a function with variables x and y"],

TextField['TF'](‘popupmenu’ = ‘PM’, ‘value’ = “cos(x^2)+sin(x*y^2)”),

["The result is"],

MathMLViewer['ML'](),

Button(“Exit”, Shutdown())

]),

PopupMenu['PM'](

Menu(“Differentiate”,

MenuItem(“With respect to x”, Evaluate(‘ML’ = ‘diff(TF, x)’)),

MenuItem(“With respect to y”, Evaluate(‘ML’ = ‘diff(TF, y)’))

),

MenuSeparator(),

MenuItem(“Exit Maplet application”, Shutdown())

)

):

> Maplets[Display](MyMaplet)

Kết quả ta được

PopupMenu

from → Maplet Programming

Soạn thảo tiếng Nga với mã UTF-8 cho MikTeX

2009 June 29
Leave a comment
tags: , , ,
by Tuan Minh

Trong phần trước có giới thiệu về gói soạn thảo tiếng Nga với mã KOI8-R. Tuy nhiên UTF-8 vẫn là mã phổ biến hơn, tiện lợi trong quá trình soạn thảo một tài liệu khoa học tiếng Nga.

Bạn cần download gói Unicode ở đây, và gói Cyrillic ở đây. Sau đó giải nén, và copy lần lượt vào các thư mục

C:\Program Files\MiKTeX 2.7\tex\latex\unicode

C:\Program Files\MiKTeX 2.7\tex\latex\cyrillic

Sau đó vào Start Menu> MikTeX 2.7> Settings> General, nhấp vào nút Refresh FNDB và Update Formats.

Để sử dụng gói này, bạn khai báo như sau:

\documentclass[a4paper,12pt]{article}

\usepackage[T2A]{fontenc}

\usepackage[utf8x]{inputenc}

\usepackage[russian]{babel}

\begin{document}

\textsc{ Я вас любил…}\\

Я вас любил: любовь еще, быть может,\\

В душе моей угасла не совсем;\\

Но пусть она вас больше не тревожит;\\

Я не хочу печалить вас ничем.\\

Я вас любил безмолвно, безнадежно,\\

То робостью, то ревностью томим;\\

Я вас любил так искренно, так нежно,\\

Как дай вам бог любимой быть другим.\\

1829\\

\end{document}

Có thể sử dụng file mẫu, download ở đây

from → TeX

« Older Entries