10-đề-tai-matlab

5/12/2018

10-đê ̀- ta i-ma tla b - slide pdf.c om

I. Phân tích ma trận theo phương pháp cholesky Thuật toán Choleski cho phép phân tích ma trận [A] thành tích hai ma trận: [A] = [L][L] T 1.Chương trình Hàm euler(): Gọi hàm : [X,Y]=euler(fxy,xo,xf,yo,n) Đầu vào : fxy:chuỗi biến với tên chức năng của tập tin chức năng có chứa f(x,y) xo,xf:giá trị ban đầu và giá trị cuối cùng của biến độc lập(vô hướng) yo:giá trị ban đầu của biến phụ thuộc tại xo(vector cột) n:số khoảng thời gian sử dụng giữa xo,xf. Đầu ra:

X = vector có chứa giá trị của biến độc lập Y = ước tính biến phụ thuộc vào mỗi giá trị của biến độc lập

2.Hàm thực hiện function [X, Y] = euler(fxy,xo,xf,yo,n) % %Giai phuong trinh y?(x) = f(x,y(x)) hay y’ = f(x) if n < 2 n = 2; end h = (xf-xo)/n; X = zeros(n+1,1); M = max(size(yo));% so phuong trinh (so cot cua ma tran Y) Y = zeros(n+1,M); %dat dieu kien dau x = xo; X(1) = x; y = yo; Y(1,:) = y';

II.

http://slide pdf.c om/re a de r/full/10-de -ta i-ma tla b

Xây dựng hàm bằng phương pháp Gauss-Jordan

1/19

5/12/2018


Mô tả phương pháp. a. Ý tưởng. Chúng ta biết rằng các nghiệm của hệ không đổi nếu ta thay một hàng bằng tổ hợp tuyến tính các hàng khác. Ta xét một hệ phương trình đại số tuyến tính có ma trận [A] không suy biến với m=n=3. Phương trình có dạng: a11x1+a12x2+a13x3=b1a21x1+a22x2+a23x3=b2a31x1+a32x2+a33x3=b3

(1) Trước hết ta khử x1 ra khỏi các phương trình, ngoại trừ phương trình đầu tiên, bằng cách nhân phương trình đầu tiên với ai1/a11 (i là chỉ số hàng) và trừ đi mỗi phương trình đó: a11(0)x1+a12(0)x2+a13(0)x3=b1(0)

a22(1)x2+a23(1)x3=b2(1)

(2)

a32(1)x2+a33(1)x3=b3(1)

Trong đó: aij(0)=aij

với i=1, j=1,2,3.

bi(0)=bi

aij(1)=aij(0)-ai1(0)a11(0)a1j(0)

bi(1)=bi(0)-ai1(0)a11(0)b1(0)

với

i,j=2,3. Việc này gọi là lấy trụ tại a11 gọi là trụ. Tiếp theo ta khử x2 trong phương trình thứ 3 của (2) bằng cách lấy phương trình thứ 2 nhân với ai2(1)/a22(1) (i=3) và trừ đi phương trình thứ 3: a11(0)x1+a12(0)x2+a13(0)x3=b1(0)

a22(1)x2+a23(1)x3=b2(1)

(3)

a33(2)x3=b3(2)

Trong đó: (4)

aij(2)=aij(1)-ai2(1)a22(1)a2j(1)

bi(2)=bi(1)-ai2(1)a22(1)b2(1)

với i,j=3

Quá trình này được gọi là thuật toán khử Gauss tiến và được tổng quát hóa thành:

aijk=aijk-1-aikk-1akkk-1akjk-1


i,j =k+1,k+2,….,m.

2/19

5/12/2018


(5) bi(k)=bi(k-1)-aik(k-1)akk(k-1)bk(k-1)

i=k+1,k+2,…..,m

b.Thuật toán Để thực hiện thuật toán khử Gauss ta dùng đoạn mã lệnh: for k=1:n-1 for i=k+1:n if A(i,k)~=0 lambda=A(i,k)/A(k,k); A(i,k+1:n)=A(i,k+1:n)-lambda*A(k,k+1:n); b(i)=b(i)-lambda*b(k); end end end Sau khi có hệ phương trình dạng tam giác ta tìm nghiệm dễ dàng. Từ phương trình thứ 3 của (3) ta có:

(6a)

x3=b3(2)a33(2)

Thay vào phương trình thứ 2 ta có: (6b)

x2=b21-a231x3a221

và cuối cùng từ phương trình thứ nhất ta có:

(6c)

x1=1a11(0)b1(0)-j=23a1j(0)xj

Ta cũng có thể tống quát hóa quá trình tìm nghiệm bằng cách tính lùi và tìm nghiệm bằng:

i=m,m-1,…,1

xi=1aii(i-1)bi(i-1)-j=i+1maij(i-1)xj

và tìm nghiệm bằng đoạn mã lệnh:

(7)

for k=n:-1:1 b(k)=(b(k)-A(k,k+1:n)*b(k+1:n)/A(k,k); end


3/19

5/12/2018


Như vậy phương pháp Gauss gồm hai bước: -khử theo thuật toán Gauss -tìm nghiệm của phương trình dạng tam giác Đoạn mà lện để tráo hàng được viết trong hàm swaprows(): function v=swaprows(v,i,j) % Trao doi hang I va hang j cua ma tran v. % Cu phap: v=swaprows(v,i,j) temp=v(i,:); v(i,:)=v(j,:); v(j,:)=temp;

Ta xây dựng hàm gauss() để thực hiện thuật toán khử Gauss function x=gauss(A,B) % Kích thước của ma trận A,B là NA x NA và NA x NB. %Hàm này dùng giải hệ pt Ax=B bằng phương pháp khử Gauss NA=size(A,2); [NB1,NB)=size(B); if NB1~=NA error(‘A và B phải có kích thước tương ứng’); end N=NA+NB; AB=[A(1:NA,1:NA) B(1:NA,1:NB)]; epss=eps*ones(NA,1); for k=1:NA %Chọn trừ AB(k,k) [akx,kx]=max(abs(AB(k:NA,k))./…


4/19

5/12/2018


max(abs([AB(k:NA,k+1:NA)epss(1:NA-k+1)]’))’); if akx1 % tráo hàng khi cần swaprows(AB,k,mx); end %Khu Gauss AB(k,k+1:N)=AB(k,k+1:N)/AB(k,k); AB(k,k)=1; for m=k+1:NA AB(m,k+1:N)=AB(m,k=1:N)-AB(m,k)*AB(k,k+1:N); AB(m,k)=0; end end %Tìm nghiệm x(NA,:)=AB(NA,NA+1:N); for m=NA-1:-1:1 x(m,:)=AB(m,NA+1:N)-AB(m,m+1:NA)*x(m+1:NA,:); end III. Phương pháp lặp Gauss-seidel Mô tả phương pháp Xét hệ phương trình dạng:


5/19

5/12/2018


Biểu diễn hệ phương trình trên dưới dạng ma trận sau: [A][x]=[b] Trong đó:

Ý tưởng chung của phương pháp Gauss-Seidel là đưa hệ Ax = b về dạng : x = Bx + g. Tức là đưa hệ phương trình trên về dạng:

Tổng quát với i=1,2,…,n ta viết:


6/19

5/12/2018


Hay viết gọn hơn là: x = Bx + g Khi đó B được gọi là ma trận lặp và:

Tiếp theo ta chọn xấp xỉ ban đầu tùy ý x10, x20, …, xn0 và tất nhiên ta cố gắng lấy chúng tương ứng với x 1, x2, …, xn (càng gần càng tốt). Tiếp theo ta giả sử rằng xấp xỉ thứ k là xi(k) của nghiệm đã biết. Theo Seidel, xấp xỉ thứ (k+1) được tính theo công thức:

Quá trình tính toán sẽ dừng lại khi xấp xỉ xk của nghiệm thỏa mãn: | xik+1 – xik |< ε , với mọi i=1,2,…n và ε sai số cho trước. Điều kiện hội tụ của nghiệm: Hệ phương trình có ma trận lặp B thỏa mãn một trong 3 điều kiện sau thì quá trình lặp sẽ hội tụ:


7/19

5/12/2018


Mở chương trình Matlab ta xây dựng hàm gauss_seidel(N,a,x0,eps,maxl) theo ý tưởng của phương pháp Gauss-Sedel. Trong đó: N: là hạng của ma trận hệ phương trình. a: là ma trận hệ số của các biến. x0: là ma trận xấp xỉ ban đầu của nghiệm (tùy chọn, tuy nhiên số phần tử phải bằng N). eps: là sai số của nghiệm tính theo Gauss-Seidel so với nghiệm đúng ( tùy theo yêu cầu bài toán). maxl: là số lần lặp lớn nhất function [x,ss,n]=gauss_seidel(N,a,x0,eps,maxl) n=0; kt=0; for l=1:maxl n=n+1;% dem so lan lap tim nghiem for k=1:N sum=0; for j=1:N if j
elseif j>k, sum=sum+a(k,j)*x0(j); end; end; x(k)=(a(k,N+1)-sum)/a(k,k);%ham lap end; ss=0; for k=1:N


8/19

5/12/2018


if abs(x(k)-x0(k))>ss ss=abs(x(k)-x0(k));%tinh sai so end; end; if ss

9/19

5/12/2018


for i = 1:maxiter r = b - a*x; x = x + omega*r; if norm(r) < tol break; end end i Để giải hệ phương trình ta dùng chương trình ctrichardsoniter.m clear all, clc a = [ 10 1 1;1 10 2; 2 2 10]; b = [12 13 14]'; x = [ 0 0 0]'; maxiter = 50; tol = 1e-6; x = richardsoniter(a, b, x, maxiter, tol) V.

Phương pháp nội suy Lagrange

function l = lagrange(x,y) %Dua vao : x = [x0 x1 ... xn], y = [y0 y1 ... yn] %ket qua: l = He so cua da thuc Lagrange bac n n = length(x)-1; %bac cua da thucl l = 0; for m = 1:n+1 p = 1; for k = 1:n+1 if k ~= m p = conv(p, [1 -x(k)])/(x(m) - x(k)); end end l = l + y(m)*p; a = x(1,1):0.2:x(1,n); hamso=polyval(l,a); plot(x,y,'o',a,hamso); end


10/19

5/12/2018


VI.

Nội suy newton

Ta xây dựng hàm newton() để nội suy: function [n,DD] = newton(x,y) %Dua vao : x = [x0 x1 ... xN] % y = [y0 y1 ... yN] %Lay ra: n = he so cua da thuc Newton bac N N = length(x) - 1; DD = zeros(N + 1, N + 1); DD(1:N + 1, 1) = y'; for k = 2:N + 1 for m = 1: N + 2 - k DD(m,k) = (DD(m + 1, k - 1) - DD(m, k - 1))/(x(m + k - 1) - x(m)); end end a = DD(1, :); n = a(N+1); for k = N:-1:1 n = [n a(k)] - [0 n*x(k)]; end Cho hàm dưới dạng bảng: x -2 -1 1 2 4 y -6 0 0 6 60 Ta dùng chương trình ctnewton.m để nội suy: clear all, clc x = [-2 -1 1 2 4]; y = [-6 0 0 6 60]; a = newton(x, y) yx = polyval(a, 2.5) VII. Phương pháp lặp đơn Giả sử phương trình (1) được đưa về dạng tương đương: x = g(x) (2) từ giá trị xo nào đó gọi là giá trị lặp đầu tiên ta lập dãy xấp xỉ bằng công thức: xn = g(xn-1) (3) với n = 1,2,....


11/19

5/12/2018


Hàm g(x) được gọi là hàm lặp. Nếu dãy xn → α khi n →∝ thì ta nói phép lặp (3) hội tụ. Ta có định lí: Xét phương pháp lặp (3), giả sử: - [a, b] là khoảng chứa nghiệm α của phương trình (1) tức là của (2) - mọi xn tính theo (3) đều thuộc [a, b] - g(x) có đạo hàm thoả mãn : g(x) <=q<1 a

12/19

5/12/2018


else fprintf('Hoi tu sau %d lan lap\n',k) end Để tính lại ví dụ trên ta dùng chương trình ctsimpiter4_2.m clear all, clc f = inline('-0.5*((x - 1).^2 - 3)'); [x, ss, xx] = simpiter(f, 0.5,.00001,200) VIII. Tính đạo hàm bậc cao Ta xây dựng hàm diffn() để tính đạo hàm tới bậc 5: function df = diffn(f, n, x) % Tinh dao ham cap n cua f tai x if n>5 error('Ham chi tinh duoc dao ham den bac 5'); return; end; N = 5; xo = x; T(1) = feval(f,xo); h = 0.005; tmp = 1; for i = 1:N tmp = tmp*h; c = difapx(i,[-i i]); %he so cua dao ham dix = c*feval(f,xo + [-i:i]*h)'; T(i+1) = dix/tmp; %dao ham end df = T(n+1); h = 0.005; tmp for i == 1; 1:N tmp = tmp*h; c = difapx(i,[-i i]); %he so cua dao ham dix = c*feval(f,xo + [-i:i]*h)'; %/h^i; %dao ham T(i+1) = dix/tmp; %he so cua chuoi Taylor end


13/19

5/12/2018


df = T(n+1); Để tính đạo hàm của hàm ta dùng chương trình ctdiffn.m clear all, clc f = inline('x.^2 + atan(x)','x'); df = diffn(f, 5, 0)

IX.

Phương pháp Euler

A. Mô tả phương pháp 1.Định nghĩa Phương pháp chuỗi Taylor cấp p=1 được gọi là phương pháp Euler xi+1 = xi + h, yi+1 = yi + hf ( x ,i yi) Về phương diện hình học,nghiệm gần đúng cho bởi phương pháp Euler được minh họa bởi một đường gấp khúc mà đoạn đầu tiên trùng với tiếp tuyến với đường cong nghiệm tại x0. 2.Phương pháp Giả sử ta có phương trình vi phân: y'x=fx,yya=α

Và cần tìm nghiệm của nó.Ta chia đoạn[xo,x] thành n phần bởi các điểm chia x0 < x1 < x2 <...< xn = x Theo công thức khai triển Taylor một hàm lân cận xi ta có:

(1) Nếu (xi+1-xi) khá bé thì ta có thể bỏ qua các số hạng (xi+1-xi)2 và các số hạng bậc cao y(xi+1) = y(xi) + (xi+1‐ xi)y′(xi)


14/19

5/12/2018


Trường hợp các mốc cách đều : (xi-1 ‐ xi) = h = (x ‐ x0)/ n thì ta nhận được công thức Euler đơn giản: yi+1=yi+hf(xi,yi) (2) Về mặt hình học ta thấy (1) cho kết quả càng chính xác nếu bước h càng nhỏ . B. Chương trình cài đặt 1.Chương trình Hàm euler(): Gọi hàm : [X,Y]=euler(fxy,xo,xf,yo,n) Đầu vào : fxy:chuỗi biến với tên chức năng của tập tin chức năng có chứa f(x,y) xo,xf:giá trị ban đầu và giá trị cuối cùng của biến độc lập(vô hướng) yo:giá trị ban đầu của biến phụ thuộc tại xo(vector cột) n:số khoảng thời gian sử dụng giữa xo,xf. Đầu ra:

X = vector có chứa giá trị của biến độc lập Y = ước tính biến phụ thuộc vào mỗi giá trị của biến độc lập

2.Hàm thực hiện function [X, Y] = euler(fxy,xo,xf,yo,n)


15/19

5/12/2018


% %Giai phuong trinh y?(x) = f(x,y(x)) hay y’ = f(x) if n < 2 n = 2; end h = (xf-xo)/n; X = zeros(n+1,1); M = max(size(yo));% so phuong trinh (so cot cua ma tran Y) Y = zeros(n+1,M); %dat dieu kien dau x = xo; X(1) = x; y = yo; Y(1,:) = y'; for i = 1:n if nargin(fxy)>1; k1 = h*feval(fxy,x,y); else k1 = h*feval(fxy,x); end y = y + k1; x = x + h; X(i+1) = x; Y(i+1,:) = y'; end function dy = f1(t,y) dy = zeros(3,1); dy(1) = y(2)*y(3); dy(2) = -y(1)*y(3); dy(3) = -0.51*y(1)*y(2); X.

Phương pháp newton(BT tối ưu)

Việc tìm điểm cực tiểu của hàm f(x) tương đương với việc xác định x để gradient g(x) của hàm f(x) bằng zero. Nghiệm của g(x) = 0 có thể tìm được bằng c dùng phươngcủa pháp Newton hệ =phương tìm nghiệm phương trìnhcho g(x) 0 là: trình phi tuyến. Hàm newtons(x) dùng function [x, fx, xx] = newtons(f, x0, tolx, maxiter) h = 1e-4; tolfun = eps; EPS = 1e-6;


16/19

5/12/2018


fx = feval(f, x0); nf = length(fx); nx = length(x0); if nf ~= nx error('Kich thuoc cua g va x0 khong tuong thich!'); end if nargin < 4 maxiter = 100; end if nargin < 3 tolx = EPS; end xx(1, :) = x0(:).'; for k = 1: maxiter dx = -jacob(f, xx(k, :), h)\fx(:);%-[dfdx]ˆ-1*fx xx(k + 1, :) = xx(k, :) + dx.'; fx = feval(f, xx(k + 1, :)); fxn = norm(fx); if fxn < tolfun | norm(dx) < tolx break; end end x = xx(k + 1, :); if k == maxiter endfprintf('Ket qua tot nhat sau %dlan lap\n', maxiter) function g = jacob(f, x, h) %Jacobian cua f(x) if nargin < 3 h = 1e-4; end h2 = 2*h; n = length(x); x = x(:).'; I = eye(n); forg(:, n =n)1:n = (feval(f, x + I(n, :)*h) - feval(f, x - I(n,:)*h))'/h2; end Để tìm cực tiểu của hàm bằng phương pháp Newtons ta dùng chương trình ctnewtons.m:


17/19

5/12/2018


clear all, clc f = inline('x(1).^2 - x(1)*x(2) - 4*x(1) + x(2).^2 - x(2)'); g = inline('[(2*x(1) - x(2) -4) ( 2*x(2) - x(1) - 1)]'); x0 = [0.1 0.1]; tolx = 1e-4; maxiter = 100; [xo, fo] = newtons(g, x0, tolx, maxiter)


18/19

5/12/2018



19/19

10-đề-tai-matlab

Recommend Documents