🔙 Quay lại trang tải sách pdf ebook Giáo trình phương pháp tính và tin học chuyên ngành Ebooks Nhóm Zalo ƠNG (Chủ biên) - NGUYỄN TRỌNG DŨNG GIÁO TRÌNH PHƯƠNG PHÁP TÍNH VÀ TIN HỌC CHUYÊN NGÀNH NGUYỄN CHÍNH CƯƠNG (Chủ biên) NGUYỄN TRỌNG DŨNG Giáo trình PHƯƠNG PHÁP TÍNH VÀ TIN HỌC CHUYÊN NGÀNH NHÀ XUẤT BẢN ĐẠI HỌC s ư PHẠM MỤC LỤC Trang Lời mỏ đầu.......................................................................................................................... 5 PHẦN I. CÁC PHƯƠNG PHÁP TÍNH VÀ MÔ PHỎNG..........................................................7 Chương 1. Các phương pháp tính và mô phỏng...........................................................7 §1. Sai số.......................................................................................................................... 7 §2. Phương pháp giải gần đúng phương trình và hệ phương trình...............................11 §3. Phương pháp nội suy hàm số ..................................................................................15 §4. Phương pháp tính gần đúng đạo hàm. Tính tích phân - Tính vi phân.................... 20 §5. Một sô' phương pháp gần đúng trong vật lí lượng tử.............................................. 26 §6. Phương pháp tính số và mô phỏng......................................................................... 37 PHẦN II. NGÔN NGỮ LẬP TRÌNH MATLAB VÀ FORTRAN............................................ 41 Chưong 2. Giới thiệu về Matlab........................................................................................41 §1. Mở đầu.......................................................................................................................41 §2. Cài đặt và khởi động Matlab 7.0............................................................................. 42 §3. Quân lí không gian làm việc của Matlab.................................................................43 §4. Một số lưu ý khi làm việc với Matlab....................................................................... 46 Chưong 3. Ngôn ngữ lập trình tính toán trên Matlab.................................................. 49 §1. Các phép toán cơ bản.............................................................................................. 49 §2. Các phương pháp gần đúng....................................................................................65 §3. M-file..........................................................................................................................76 §4. Symbolic Math Toolboxe...........................................................................................81 ChƯdng 4. ĐÕ hoạ............................................................................................................ 110 §1. Đổ thị hai chiều....................................................................................................... 110 §2. Đồ thị ba chiều........................................................................................................116 §3, Một vài dạng đồ thị đặc biệt...................................................................................119 3 Chương 5. Mô phỏng Simulink Power System Blockset và G ui............................126 §1. Mô phỏng bằng Simulink...................................................................................... 126 §2. Mô phỏng bằng Power system blockset............................................................. §3. Gui trong Matlab....................................................................................................148 §4. Thiết kế phần mềm mô phỏng trên Matlab.......................................................... 167 Chương 6. Một sô' ứng dụng của ngôn ngữ lập trinh Matlab..................................184 §1.Môphỏng một số cơ hệ.........................................................................................184 §2. Mô phỏng một số bài toán phần điện - từ............................................................ 202 Chường 7. Ngôn ngữ lập trinh Fortran 9.0................................................................. 211 §1. Cấu trúc chương trình Fortran............................................................................... 211 §2. Các cấu trúc lệnh lặp của Fortran.........................................................................218 §3. Chương trình con....................................................................................................228 §4. Mảng.......................................................................................................................237 §5. Xây dựng cơ sà dữ liệu.......................................................................................... 240 §6. Kiểu file...................................................................................................................247 §7. Một số ví dụ ứng dụng của Fortran trong vật l i ....................................................256 Tài liệu tham khảo...............................................................................................................263 4 LỜI M ỏ ĐẨU Trong quá trình nghiên cứu, học tập các môn học chuyên ngành Vật lí, chúng ta thường gặp phải một sô' khó khăn khi tính toán như: Tính một tích phân khỏng ở dạng thông Ihường, giải một phương trình vi phân nhiều ẩn, biến đổi một biểu thức vật lí có dạng toán học phức tạp, nghiên cứu một hệ vật lí giả định... Việc này rất tốn công sức và không phải lúc nào cũng có thê thu được những kết quả tường minh. Vì vậy, đòi hỏi chúng ta phải sử dụng phương pháp tính số, mô phỏng,... từ đó có thể tìm lại các kết quả ờ dạng gẩn đúng với sai số nằm trong giới hạn cho phép. Để làm được điểu đó thì ngoài những kiến thức chuyên ngành tôì đòi hỏi chúng ta phải có kiến thức về phương pháp tính số và các lĩnh vực khoa học tính toán mô phòng. Cùng với sự phát triển nhanh chóng cùa ngành Công nghệ thòng tin, việc sử dụng máy tính và các phần mềm chuyên dụng đã và dang là xu thế tất yếu cùa các ngành khoa học tính toán nói chung và ngành Vật lí học nói riêng. Thực tế chì ra rằng, việc sù dụng máy tính điện tử trong nghiên cứu, học tập vật lí đã đem lại những thành tựu vô cùng quan trọng. Điều này được thể hiện ờ việc ứng dụng máy tính trong rất nhiểu chuyên ngành, nhiều ĩinh vực khác nhau của vật lí như: Vật lí lí thuyết, Thiên văn, Vật lí chất rắn... nhằm hỗ trợ cho việc tính toán và đo đạc í chính xác, xây dựng các mô hình lí thuyết, mô phỏng, dự đoán các hiện tượng và hỗ trợ thí nghiệm... Việc sử dụng máy tính nhằm hỗ trợ quá trình nghiên cứu đã trờ nên phổ biến đối với công tác nghiên cứu khoa học trong nước và trên thê giới. Từ mong muốn đó, giáo trình 'Thương pháp tính và tin học chuyên ngành" bước đầu trang bị cho học viên cao học ngành Vật lí phương pháp tính gần đúng, tính số và ngôn ngữ lập trình Mathlab, Fortran... nhằm hỗ trợ mô phòng các hệ vật lí. Đây là những công cụ được nhiều nhà khoa học, học viên, sinh viên một số trường đại học tiên tiến sử dụng. Phương pháp tính số và mô phỏng là một lĩnh vực rộng và mới mẻ nên có nhiều vấn đề cần quan tâm, do thời gian còn hạn ch ế nên trong giáo trình khó tránh khỏi những thiếu sót. Rất mon® được quý bạn đọc góp ý để chúng tòi chinh sửa và hoàn thiện giúp cho các học viên cao học và các sinh viên có được một tài liệu học tập hữii ích. CÁC TÁC GIẢ 5 P H Ầ N I CÁC PHƯƠNG PHÁP TÍNH VÀ MÔ PHỎNG Chương 1 CÁC PHƯƠNG PHÁP TINH VÀ MỒ PHỎNG §1. SAI SỐ Sai s ố thường xuyên xảy l a khi chúng ta chuyển bủi toán phức rạp về dạng đơii giàu hơn. Mặt khác, sai sô'l òn phản ánli kết quả tính toán có sát với các giú trị thực liuy klìông. 1.1. Sai số tuyệt đối và sai số tưdng dối Sai số tuyệt đôi: Xét đại lượng đúng A có giá trị gần đúng là u. Lúc đó ta nói “a xấp xỉ A” và viết a « A . Trị tuyệt đối |a - A| gọi là sai số tuyệt đối của a. Do không biết số đúng A nên không tính được sai số tuyệt đối của a, vì vậy ta tìm cách ước lượng sai số đó bằng số dương A„ nào đó lớn hơn hoặc bằng ịa - Aị a - A < Aa ( l . I) Số dương An này gọi là sai số tuyệt đối giới hạn cùa a. • Nếu A„ là sai số tuyệt đối giới hạn cúa a thì mọi số A'> Aa đều là sai số tuyệt đối giới hạn cùa a. Vì vậy trong điều kiện cụ thể người ta chọn A.J là số dương nhỏ nhất có thể được thoà mãn ( I. I). • Nếu số xấp xỉ a cùa A có sai số tuyệt đối giói hạn là A(, thì quy ước viết: A = a ± A., Với nghĩa của ( I. I) tức là: Ỉ1 — A., ^ A ^ ÍI + A., ( 1.2) (1.3) Sai sốtưưng đôi: Ti số|a - Aj |a - agọi là sai số tương đối cùa a so với A. la A 7 Ta gọi ti số: (1.4) a là ti sô tương đối giới hạn của a. Từ đó ta có : A a = | 4 5 a Do (1.5) nên ( 1.2) cũng có thế viết: A = a .(l± õa ) (1.5) ( 1.6) * Lưu ý: Sai số tuyệt đối khỏng nói lên chất lượng cùa một số xấp xỉ, “chất lượng” ấy được phàn ánh qua sai sô tương đối. Cho nên khi muốn so sánh sự chính xác cùa các phép đo cùng loại thì ta sừ dụng sai số tương dối. 1.2. Cách viết số xấp xỉ Chữ số có nghĩa: Một số thập phân có thể gồm nhiều chữ số nhưng ta chi tính các chữ sô' khác không đầu tiên từ trái sang phải là chữ số có nghĩa. Chữ sô đáng tin: Mọi số thập phân đều có dạng: a = ± X a s.10s (I.7) trong đó: a s là những sô' từ 0 đến 9. Già sử a )à giá trị xấp xỉ của A với sai số tuyệt đối giới hạn Aa . Nếu Aa < 0.5.10s thì nói a s là chữ số đáng tin, nếu Aa> 0 .5 .l0 s thì nói a s là chữ sô'đáng nghi. • Nếu a s là đáng tin thì tất cả những chữ số có nghĩa đứng bên trái nó cũng là đáng tin. • Nếu a s là đáng nghi thì tất cả những chữ số có nghĩa ở bên phải nó cũng là đáng nghi. Cách viết số xấp xi: Cácli 1: Viết kèm theo sai sô' như công thức (1.2) hoặc (1.6). Cách 2: Viết theo quy ước “Mọi chữ sô' có nghía đểu đáng tin”. Có nghĩa là sai số tuyệt đối giới hạn khôns lớn hơn một nứa đơn vị cùa hàng cuối cù n s. 1.3. Sai số quy tròn Hiện tượng quy tròn và sai số quy tròn: Trong tính toán khi gặp một số có quá nhiều chữ số đáng nghi người ta bó di một vài chữ sò' cuối cho gọn. việc làm đó gọi là quy tròn số. Mỗi khi quy tròn một số người ta tạo ra một sai số mới oọj là sai số quy tròn, nó bằng hiệu giữa số quy tròn và số chưa quy tròn. 8 Quv tác: Quy tròn sao cho sai số quy tròn tuyệt đối không lớn hơn một nừa đơn vị ờ hàng được giữ lại cuối cùng, tức là 5 đơn vị ờ hàng bò đi đầu tiên. Cụ the là: Nếu chữ số bỏ đi dầu tiên không nhỏ hơn 5 thì thêm vào chữ số giữ lại cuối cùng một đơn vị, còn nếu chữ số bỏ đi đầu tiên nhỏ hơn 5 thì để nguyên chữ số giữ lại cuối cùng. Sai số đã quy tròn: ( 1.8 ) 1.4. Sai số tính toán và sai sô' phương pháp Khi giải gần đúng một bài toán phức tạp ta phải thay bài toán đã cho bằng một bài toán đơn giản hơn có thể giải được thông qua việc thực hiện các phép tính thông thường bằng tay hoặc trên máy tính. Phương pháp thay bài toán phức tạp bằng bài toán đơn giản như thế gọi là phương pháp gần đúng. Sai số do phương pháp gần đúng tạo ra gọi là sai sô'phương pháp. Khi giải các bài toán đơn giản, và thực hiện các phép tính thòng thường ta luôn phải quy tròn các k ít quà trung gian. Sai số tạo ra bời tất cả các lần quy tròn như vậy gọi là sai số tính toán. * Sai số cuối cùng là tổng hợp cùa hai loại sai số phương pháp và sai số tính toán. 1.5. Các quy tắc tính sai số Xét hàm số u của hai biến số X và y: u = f(x,y) Cho biết sai số về X và y, ta cần lập công thức tính sai số về u. Kí hiệu: Ax, Ay, Au chỉ các sô' gia cùa X, y, u. dx, dy, du chi’ các vi phân của X, y, u. Ax , Ay , Au là sai sô' tuyệt đối của X, y, u. Ta luòn có: |Ax| < Ax; |Ay| b n>... tạo nên dãy đơn điệu không tăng và bị chặn dưới bởi số a, nên khi n -> +00 tìr (1.25) ta nhận được: lim an = lim bn = a n—>+«: 11—»+« Cho II -» +00 trong bất đẳng thức (1.24), do sự liên tục cùa hàm số f(x) ta có: |f(a )|' < 0 Vậy f(a ) = 0 và a là nghiệm của phương trình. Đánh giá sai số của nghiệm gần đúng: Trong thực hành ta không Ihế thực hiện phương pháp chia dôi vô hạn lẩn để nhặn được nghiệm đúna của phương trình mà chi có thể áp dụng n lẩn phương pháp chia đòi, với n là một số nguyên dưona hữu hạn. Dừng lại ờ lần thứ 11 ta có : a n < c c < b n và b n - a n = — ( b - a ) 2 " 13 Phương pháp này có ưu điểm là đơn giản, dễ lập chương trình chạy máy tính nhưns lại có nhược điểm là tốc độ hội tụ chậm. 2.3. Phương pháp lặp Mô tả phương pháp: Giả sử phương trình f(x) = 0 có nghiệm thực a phân li trong khoảng [a, b]. Trước hết ta chuyên phương trình f(x) = 0 về dạng: x = cp(x) (1.26) Sau đó ta chọn một số X g nào đó 6 |a, b] làm xấp xỉ đầu rồi tính dần day sô' .v„ theo quy tắc: x n = < p (x„-i) n = 1, 2,... (1.27) Quá trình tính này lặp đi lặp lại nên phương pháp ờ đày là phương pháp lặp và ọ(x) là hàm lặp. Sự hội tụ của phương pháp lặp: Nếu dãy x n —> a khi n —>+00 thì ta nói phương pháp lặp (1.26), (1.27) hội tụ. Định lí 3: Xét phương pháp lặp (1.26), (1.27) giá sừ: 1) [a, b] là khoảng phân li nghiệm a cùa phương trình f(x) = 0 2) Mọi x n tính theo (1.26), (1.27) đều e [a, b] 3) Hàm ọ(x) có đạo hàm thoả mãn |cp'(x)Ị < q < 1, a < X < b (1-28) trong đó q là một hằng số. Thế thi phương pháp lặp (1.26), (1.27) hội tụ. * Nếu (p’(x )> 0 ta có thể chọn XH 6 |a. b] một cách bất kì, nếu (p'(x) < 0 thì phái chọn X0 theo quy tắc sau: 14 , a + b , x0= b k h i —-— = c < a < b (1.29) Muốn biết a thuộc nửa khoảng nào ta chỉ việc tính f(c) rồi so sánh dấu cùa nó với dấu của f(a) Đánh giá sai số + Công thức đánh giá sai số thứ nhất | a - x n| < - 9 - | x n ( 0 m > 0 . BÀI TẬP (1.31) 1.4. Giải gần đúng phương trình: sinx + X2 = 1 bằng phương pháp chia đôi với kết quả có hai chữ số đáng tin. 1.5. Giải gần đúng phương trình: cosx - X = -0,15 bằng phương pháp chia đôi với kết quả có ba chữ số đáng tin. §3. PHƯƠNG PHÁP NỘI SUY HÀM s ố 3.1. Đa thứ c nội suy Vàn đề nội suy: Trong thực tế nhiều khi phải phục hồi một hàm số f(x) tại mọi giá trị cùa X trẽn đoạn [a, b] mà chỉ biết một số hữu hạn giá trị cùa hàm số tại một số hữu hạn các điểm rời rạc cùa đoạn đó. Các giá trị đó được cung cấp qua thực nghiệm hay tính toán. Vậy nảy sinh một vấn đề toán học sau: Trên đoạn |a, b| cho một lưới các điểm chia Xj, i = 0, l n : a < X „ ,X |,...,x n < b và tại các nút X| cho giá trị của hàm số y = f(x) là: 15 y, = f(X ị), i = 0, 1,..., n Ta viết thành báng sau: X *0 X| 1 y yo y. y„ (1.32) Hãy xây dựng một đa thức bậc n: pn(x) = a„x" + ii|.xn +... + an_| .X + a n Sao cho p n(Xị) = y j, i = 0, 1,.... 11. Đa thức p n(x) gọi là đa thức nội suy cùa hàm f(x). Ta chọn đa thức nội suy hàm f(x) vì đa thức là loại hàm đơn gián, luôn có đạo hàm và nguyên hàm, việc tính giá trị cũng đơn gián. Ta có: pn(x) = (...(a().x + a, )x + a 2).~) + a n bo = ao’ b| = by.c + U|, Sơ đồ Hoocne lính giá trị p n(c): Sự duy nhất của đa thức nội suy: bi = b| .c + a->. (1.33) b n = b n - l - c + a n = P n ( c ) Địnli l í 5: Đa thức nội suy pn(x) cùa hàm sô f(x) định nghĩa như trên nếu có thì chí có một mà thôi. * Như vậy ta thấy đa thức nội suy có thê xây dựng nhiều cách, nhimg vì nó có tính duy nhất nên tất cá các dạng của nó đều có thê quy về nhau được. 3.2. Đa thức nội suy Lagrange Thành lập đa thức nội suy Lasrange: , (x) = (x -x „)....(x -x ,_ , ) ( x - x it,)...(x -x „ ) (Xị - x „ ) . . . ( x i — Xị_, )(Xị - x ^ , ) . . . ( x : - x „ ) Rõ làng /, (X) là đa thức bậc n và ta có : (1.34) (1.35) 16 1 khi j = i 0 khi j * i Ta gọi /,(xj) là đa thức Lagrange cơ bán. Bày giờ xét biểu thức: p,(x) = I y / ( x ) Ta thấy p„(x) vừa là một đa thức bậc n, vừa thoả mãn: p n(Xị) - y i. Ta gọi nó là đa thức nội suy Lagrange. Nội suy bậc nhất: Với n = 1 đa thức nội suy sẽ là: fp,(x) = y0/0(x) + y,/,(x) • , , , x - x , x - x „ O -3 6 ) X, - x 0 Do đó: P |(x ) = y 0 . - — + y , . x - x „ Xc - X ị x , - x 0 Đa thức P | (x) là bậc nhất với X có dạng: Av + B Nội suy bậc hai: Với n = 2 đa thức nội suy có dạng là: P i(x) = yI)-/0(x) + y r/1(x) + y,./,(x) (1.37) Trong đó : /0(x) = /,(x) = /2(x ) = ( x - x ,) .( x - x 2) (x0-x ,) .(x 0- x 2) ( x - x 0) .( x - x 2) (X, - x 0) .( x ,- X , ) (x - x „ ) .( x - x ,) (x2- x 0).( x ,-x ,) Đa thức p->(x) là một đa thức bậc hai đối với X có dạng Ax: + Bx + c. Sai số nội suy: Định lí 6: Nếu hàm f(x) liên tục trên [a, b] và có trong (a, b) đạo hàm liên tục đến cấp n + 1 thì sai số nội suy có biểu thức : rn(x) = f (n+1)(c) n (x ) ,c e [a ,b ] (n + 1)! trong đó: n ( x ) = (X - x 0 ).(x - X, )....(x - x n ) 3.3. Đa thức Niutơn Đa thức nội suy Niutơn tổng quát: + Xây dựng đa thức nội suy Niutơn xuất phát từ nút x„ Giá sử hàm V = t(\) có bàns giá trị: X Xo X| \„ y y» y> y„ (1.38) (1.39) 17 Tỉ hiệu cấp một của y tại Xj,Xj là: y [ x j , x j ] =y i- Y i Xj - Xj Tỉ hiệu cấp hai cùa y tại x j , x j , x k là: y£xị,xj,xkj Các tỉ hiệu có tính chất đối xứng: y [x i,x j] = y [x j,x ,] và y [x i,x J,x k] = y [x k,x j,x i] ( X i - X k ) Với y(x) = p„(x) là m ội đa thức bậc n thì ti hiệu cấp một tại X, x0 là : [x, J . P .M - P .O O P n [x ,x o] = ( x - x 0) là một đa thức bậc n - 1. Tỉ hiệu cấp hai tại X, x„, X, là: rc bậc n - 1. Tỉ hiệu cấp hai tại X, x„, , fr r ir 1 P n k x 0] - p n[x0,x ,] D [ X , X0 ,X |J - ;( X - X |) La---o 'TÍA'., r*tí là một đa thức bậc n - 2. Tiếp tục xây dựng tới cấp n + 1 thì: pn[x,xo,x ,,...,x n]= 0 Từ định nghĩa tì hiệu ta suy ra : p n(x ) = p n(x 0 ) + ( x - x u).pn[x ,x 0 ] Pn[X’Xo]= P n[xO’Xl] + (X- x l)Pn[x>Xo’Xl] p n[x ,x o , x | ] = p n[xo, x 1, x 2] + ( x - x 2 ).p n[x ,x o , x | , x 2 ] p n[x,x0, x ll...,x n_1] = p n [x0,x 1,x 2,...,x n ] + ( x - x n ).pn [x,x0, x |, x 2.....x n] Vì pn[x,x0,x |,...,x n] = 0 , nên từ đó ta có: p „ ( x) = P „ ( x J + ( x - x J . p „ [ x 0,x,] + (x-x 0) . ( x - x , ) . p n [x ,x , x 1, x 3] + ... (1.40) + (x - x „ ) ,( x - X, ) ....(x - x n_ j ) p n [ x ,x X ] Nếu p„(x) là đa thức nội suy của hàm y = f(x) thì: p„(x,) = f(x,). i = 0, 1.....n. Và ta có: pn(x) = y (x o) + ( x - x 0).y [x 0, x ll + ( x - x o) . ( x - x l ).y fx ,x „ ,x .,x J + ... + ( x - x o) ,( x - x l)....(x-x„_l)y[x,xll.....X„J(1.41) IX Đa thức này gọi là đa thức N iu tm ứến xuất phát từ nút xc cùa hàm y = f(x). + Xây cliữig đa thức nội suy Niutonì xuất pliúr từ I l l 'l l x„: Tương tự cách xây dựng như Irên ta có: p„(x) = y(x„) + ( x - x n).y[x„>x „-i] + (X - x n ).(x - x n_, ).y [x,x „ , x„_,, x„_, ] + ... (1.42) + ( x - x n) , ( x - x 0_1) . . . . ( x - x , ) y [ x , x n,. .. , x0] Đây là đu thức Niutơi! lùi xuất phát từ nút x„ của hàm y = f(x) Trường hợp cấc nút cách đều Xị = x0 + i.h, i = 0, n. Khi đó h gọi là bước nội suy. + Khái niệm sai phân tiến: Sai phân tiến cấp một tại i: Ay ị = y i+| - y ị Sai phân tiến cấp hai tại i: A2ys = A(Ayị) = y i+, - 2 . y +1+Jíị Sai phân tiến cấp n tại i là: Anyj=A (A n 'y ,) Khi đó ta có : y ị x ^ x , ^ - ^ 2- y[x0,x ,,x ,] = y ^ y[*0.< ......* „ ] - % n!h Bây giờ đặt X = xc + h.t trong đa thức Niutơn tiến (1.41) ta được : . t(t — 1) .2 t.(t - l)...(t - n + 1 ) n ,, P n W L = x 0 + h, = y 0 + t . A y 0 + - ^ p L A y» + ••■ + ------------ ^ ------------- A y 0 (1.43) gọi là đa thức Niutơn tiến xuất phát từ x0 trong trường hợp các nút cách đều. + Khái niệm sai phân lùi Hoàn toàn tương tư như phân sai tiến ta có: Ayi = y j- y ,- i Any i =A(An- ,y i) 19 Ta có đa thức nội suy Niutơn lùi xuất phát từ X trong trường hợp nút cách dểu: Pn(x)lx=x0+h, = y» + t.Ay0 + A2y0 + ... + H i ~ I M . - n t i i Anyu (1-44) BÀI TẬP 1 1.6. Tính tích phân: Jsin (x 2).dx bằng cách lập đa thức nội suy gồm 6 điểm cách 0 đều tính từ điểm X = 0 của hàm sin(x2) rồi tính tích phân. 1.7. Một vật khối lượng m = 10"" (kg) được gắn vào con lắc lò xo có độ cứng k = 200 N/m. Đầu kia của con lắc được gắn chặt vào tường, hệ dao động dọc theo trục Ox trên một mặt phẳng nhẵn nằm ngang với hệ sô ma sát Ĩ1 = 0,05. Hệ tọa độ dược chọn sao cho vật m khi ờ trạng thái cân bằng nằm trùng với gốc tọa độ. Tại thời điểm ban đầu người ta cung cấp một động năng cho vật, sau đó đo đạc thực nghiệm xác định vị trí của vật có kết quả theo bảng sau: t (s) 0 0,3 0,7 0,9 1,2 X (cm) 0 0,7 1,95 1,8 0,3 Xác định vị trí, vận tốc, gia tốc của vật tại thời điểm t = ls. §4. PHƯƠNG PHÁP TÍNH GẦN ĐÚNG ĐẠO HÀM TÍNH TÍCH PHÂN - TÍNH VI PHÂN 4.1. Tính gần đúng đạo hàm Áp dụng đa thức nội suy: Đế tính gần đúng đạo hàm của hàm f(x) tại X tức là f’(x) ta có thể thay hàm f(x) bằng đa thức nội suy p(x) rồi tính đạo hàm cùa đa thức nội suy: p '(x ), lấy p'(x) là giá trị gần đúng của f ’(x). Cách tính này chí có ý nghĩa khi các |xi+| — XjI nhỏ vì nếu khôna thì sai sô có thể rất lớn. Áp tlụng công thức Ta.vlor: Theo côns thức Taylor ta có: f(x + h) = f(x) + h.f'(x) + — f ”(c) (c = X + e.h, 0 < 0 < I ) 20 Khi |h| bé thì số hạng cuối ở vế phải rất bé, ta có thể bỏ qua và có: f(x + h ) - f ( x ) » h.f'(x) (1.45) f(x + h ) - f ( x ) Vậy ta có: f'(x )- (1.46) 4.2. Tính gần đúng tích phân xác định b Xét tích phân xác định: I = Jf(x)dx a Nếu f(x) liên tục trên [a, b] và có nguyên hàm là F(x) thì công thức Niutơn - Lépnit cho: b Jf(x)dx = F (b )-F (a ) a Nhưng nếu không tìm được nguyên hàm của f(x) ở dạng sơ cấp hoặc nguyên hàm đó quá phức tạp thì tích phân I phải tính gần đúng. Sau đây ta sẽ trình bày hai công thức tính gần đúng tích phân I dựa trên tư tường thay hàm f(x) bằng một đa thức nội suy. Công thức hình thang: Chia [a, b] thành n đoạn cong bằng nhau bời các điểm chia Xji a = x0< X, <... < xn_, < x„ = b (Xị= a + i.h , h = ^ ~ a\ i = 0, 1,..., n). Đặt y { = f(x,) ta có: bb X, x 2 n Ịf(x )d x = j f d x + Jfd x + ...+ jfd x a x0 X | n-1 (1.47) Để tính tích phân ờ vế phải ta thay hàm f(x) bằng một đa thức nội suy bậc nhất p,(x). Với tích phân thứ nhất ta có: Jf(x )d x = Jpi (x)dx x0 x0 Đổi biến X = x0 + h.t, ứng vởi x0 là t = 0, ứng với X, là t = 1 nên có: 1=1 *1 I t 2 Jp i (x)dx = h. J(y 0 + t.Ay 0 )dt = h.(y,)t + ■— Ay0) , Ayo = h. y « + - y L x 0 0 _ h y ọ + y i 1=0 21 Vậy có: X, X, = h yọ + yi |f(x )d x = jp ;(x )d x = h 2 »0 *0 Về mặt hình học có điều đó có nghĩa: thay diện tích hình thang cong XqMqMiX, bời diện tích hình thang thường X0M0M,X| Đối với tích phân thứ i + 1 ta có: y \ + Y i + I jf(x)dx a h,- Vậy công Ihức (1.47) cho: b Jf(x)dx = ệ .[(y 0 + y 1) + (y l + y 2) + - + (y n-i + y „))] (y<> + yn)+ yi + y 2 + —+ y n-i trong đó: h = b - a Cônq thức Iiàv ?ỌÍ là công thức hình thunẹ. Đánh giá sai số: Người ta chứng minh được: | I - I T|< — ,h2( b - a ) I TI 12 M = max f (x), a < X < b (1.48) (1.49) Công thức Simxơn: Ta chia [a, b] thành 2n đoạn cong bằng nhau bời các điém chia x:: a = x„ < X, <...< x;„ = b _ g Xi = a + i.h. h = ------ , i = 0. 1,.... 2n. Giả sử y, = f(x,). Ta có: b \ _< \ V, jf(.\)dx = JVd\ + JYclx+...+ I fclx (1.50) li Để tính mỗi tích phân ờ vế phải ta thay f(x) bằng đa thức nội suy bậc hai. Với X-> X-> tích phân thứ nhất ta có: Jf(x)dx = Jpj(x)dx *0 *0 Đổi biến X = x0 + h.t thì dx = h.dt, ứng với x0là t = 0, úng vói X, là t = 2. Do đó: x 2 2 J p 2(x)dx = h .J(y 0 + tA y0 + ^ ^ A 2y 0 )dt = - ( y 0 + 4 yi + y 2) *0 ố Vậy ta có: *2 h J f ( x ) d x = j( y 0 + 4 y 1 + y 2) Xo Đối với tích phân sau cũng có cách tính tương tự: J f(x)dx = | ( y 2i + 4 y2i+1 + y :i+2) x2i Vậy ta cũng được: 1 = ^ -í(y o + y 2 n )+ 4 (y i + y 3 + - + y 2 n - i ) + 2 ( y 2 + y 4 + - + y 2 n - 2 ) ] . (151) 3 với h = b - a 2n Còng thức này gọi là công thức Sim.xơii. Đánh giá sai số: Người ta chứng minh được: | I - I S|< — .h4( b - a ) 1 1 180 (1.52) M = m ax|f<4)(x)|, a < X < b 4.3. Phương pháp giải phương trình vi phân Bài toán Côsi đối với phương trình vi phân cấp một: Cho khoáng [x„. X|. Tìm hàm số y = y(x) xác định trên [x,„ X | và thoá mãn: y' = f(x.y). X(| < X < X (1.53) y(X()) = n (1.54) 23 trong đó fix, y) là m ột hàm sô đã biết cùa hai đối sô' X y còn rj là một số thực cho trước. PliươìiiỊ pháp clmỗi Taylor: Ta tìm nghiệm y(x) khai triển thành chuỗi Taylor tại X = x0: y(x)=y(x0) + ^ ^ . ( x - x <)) + ^ p ( x - x 0)2+....+^— (x -x 0)k +... (1.55) 1! k! Tính các đạo hàm y(k)(xu) cùa y tại x„ ta có: y’(x0) = f(x0, y(x„)) (1.56) y'(x) = (y')' = (f(x,y(x)))' = -^ (x ,y (x )) + — (x,y(x)).y'(x) (1.57) ổx ỡy Thay X = x0 và y(x0) = r| ta được: y”(x0) = (y')' = (f(x 0,ĩl))' = |^ ( x o .ĩ|) + |^(Xo,r|)f(Xa,Ti) ổx ổy Một cách tưcmg tự, để tính y”(x0): trước hết ta phải lấy đạo hàm (1.57) sau đó thay X = x0 và cứ thế tiếp tục. Với X khá gần x0 thì chuồi Taylor hội tụ về nghiệm của bài toán và độ chính xác càng cao khi n càng lớn. Phương ph áp ơ le : Chia đoạn [x0, Xị thành n đoạn nhỏ bằng nhau bời các đ iể m X; X = x„ + i.h, i = 0, 1,..., n h X - x 0 (1.58) n Phương pháp ơle cho phép tìm cách tính gần đúng giá trị cùa y(x) chỉ tai các nút X, mà thôi, chứ không phải tại m ọi x e [x(), x ] . * Cọi y(x) là nghiệm cùa bài toán Côsi và y(x,) là giá trị của y (x ) tại X LI, là giá trị gần đú ng cù a y (x ,) m à ta m u ốn tính. Sau đây ta xây dựng còng thức tính 11,. Giá sử đã biết 11, tại nút X, và muốn tính u,+, tại nút x,+l. Khai triển Taylor: y(x)=y(xi)+ -^Y ^.(x-xi) + ^ - ^ ( x - x i)2 c ,= x + 0 ( x -x ).O < 0 < 1 Thay X = x,*! = X, + h và y'(x,) = f(x,, y(x,)), ta được: 24 y(xi+|) = y(xi)+K y'(Xj) , h2 y*(ci) 1! 91(1.59) Khi h bé, số hạng ờ cuối vế phải có thể xem là bé, không đáng kể ta bỏ qua và thay y(Xj) bằng Uj ta được công thức: ui+l= u i + h i.f(xi>ui) (1.60) Công thức này cho phép tính ui+l khi biết Uị và điéu kiện Côsi đặt u0 = r |. Sự hội tụ của phưcmg pháp ơle: Ta gọi e, = Uj —y ( X j) là sai số phương pháp ơle. Nếu tại Xị xác định, e; -> 0 khi h -> 0 , tức là Uj -> y(Xj) khi h -> 0 thì ta nói phương pháp ơ le hội tụ. 4.4. Phương pháp giải hệ phường trình vi phân Hệ phương trình: Cho khoảng [x0, X]. Tìm hai hàm số y = y(x) và z = z(x) xác định trên đoạn [x0, X] và thoả mãn y ’ = f(x,y,z), z’ = g(x,y,z) y(x o) = Tli > z(*o) = 1l 2 trong đó: T|i,T|2 là hai số thực cho trước. (1.61) (1.62) Phương pháp chuỗi Taylor: Giống như đối với phương trình ta cũng có thể mờ rộng áp dụng cho bài toán hệ. Với phương pháp Taylor ta viết: y ( x ) = y ( x o ) + ^ ^ .( x - x 0) + ^ - ^ ( x - x 0)2 + ....+ ^ -Ị-(x -x 0)k +... 2(x) = z(x 0) + ^ . ( x - X 0) + ^ ( x - X 0)2 + ....+ ^ ( x - X 0)k + ... Với: y(x0) = T|1, z(x0) = t|2 y ’(x0) = f(x0, y(x0), z(x0)) = f(x 0,r|i,T l2) z ’(x0) = g(x0, y(x0), z(x0)) = g (x0,T ii,n2 ) Phương pháp ơ le: Trước hết ta chia đoạn [x0, X] thành n đoạn con, giả thiết bằng nhau cho đơn giản, bời các điểm: n Phương pháp ơle viết: Uj + | = U| + hj -f(Xj ,Uj , Vj) 25 vi+l = V, uo= n ,> vo = rb Như vậy biết u,, Vị tính được ngay u,t |. v,+|. BÀI TẬP 1.8. Cho tích phân 1 ='c dx y2+C09Í Bằng cách chia đoạn [0, 1 ] thành 10 đoạn bằng nhau, hãy tính gần đúng I theo a. Đa thức nội suy. b. Công thức hình thang. c. Công thức Simxơn. 1.9. Tính tích phàn I = je“2x'c dx 0 bằng cách khai triển hàm dưới dấu tích phân thành chuỗi luỹ thừa. §5. MỘT SỐ PHƯƠNG PHÁP GẦN ĐÚNG TRONG VẬT Lí LƯỢNG TỬ Khi ngliiên cừu hệ lượng tử có tươìig tác thì Hcimiltoiiian cùa hệ hrcnig tử buo gồm số hạng klìôiig tươiìg tác và số liạng tưcniạ tác. Các hài toán có lien quan đẽh số hạng tương tác hầu hết kliông thể giải dược một cách cliinh xác. Các sô'hạng tương tác này gọi là các nhiễu loạn gáy bởi trưởng Iiẹoài hay ilo tưcnig tác cùa chính hệ lượiig từ ẹũy ra. 5.1. Phương pháp nhiễu loạn đối ■với trạng thái dừng không suy biến Giá sử Hamiltonian của hệ lượng tử có thế viết dưới dạng: H = Hu +ÀHi ( I.63) trong đó: + Hí) là Hamiltonian cùa hệ ờ trạng thái không có nhiễu loạn, với trị riéng EỊ,'” và veclơ riêns l^ n ) khône suy biến đã được biết chính xác: 26 HoK ) = E<;>K> (1.64) + Hi là thành phần nhiễu loạn, X là thông sô' nhiễu loạn. Kí hiệu các vectơ riêng của H là |*Pn) ứng với trị riêng E n , ta có: (H„ + m 1) |4 ' n) = En|4 'n) - (1.65) Khai triển theo các vectơ trạng thái riêng |„) và giả thiết có sự tương ứng giữa các trạng thái ^Fn) và I ). Nghĩa là khi Ằ. 0 ta có sự tương ứng sau: |T n) ^ |4 > n) và E n -» E (0) I ' O H ' O + I c n k W K ) k* n ( 1.66) Từ điều kiện X —> 0 ta có I'-Fn) —>• |O n)dẫn đến cnk(0) = 0 . Các hệ số khai triển cnk và trị riêng E n cũng được khai triển chuỗi theo Ả. cnk( \) = x l " t + x2.c % + x \ c % + ... En = E(n0) +X.EÍ,0 +A.2.E(n2) +A.3.E(„3) +... (1.67) Thay các khai triển vào phương trình Schrodinger ta được: (H o + ú ìi) = ( E(0) +X..E(nI> +X.2.E’2> (1.68) Đây là phương trình thoả mãn với X bất kì, các hệ số mỗi bậc cùa X ở cà hai v í của phương trình phải bằng nhau. Trong gần đúng bậc một, chi giữ lại các số hạng ti lệ với X, phương trình ( 1.68) có dạng: H o X c"> 10k) + H, I „ I và sử dụng điều kiện trực giao của các vectơ trạng thái không nhiễu loạn („ |k M n k ta tìm được phần bổ chính bậc một của năng lượng: E11’ =(„ |H ,|® „) (1.71) Tương tự, nhân trái hai vế cùa phưcmg trình (1.70) với (O m I (m ^ n ) và sử dụng điểu kiện trực giao của vectơ trạng thái không nhiễu loạn ta dược: + (E^> - E'n0) ).<■„*, = 0 Ta suy ra hệ số khai triển trong gần đúng bậc một: > K Ị h . K ) nk do đó trạng thái trong gần đúng bậc một có dạng: k,„ E„ - Ek Bỏ qua các số hạng chứa luỹ thừa của bậc ba trờ lên, ta dược: (1.72) (1.73) Ho ciỉ' I ■k) + H, c<;> I = E<:> c«' I IOk) + E<„" cl1 > I o k) + E‘„2) I „) (1.74) k*n Nhân trái vô hướng hai vế của (1.74) với bra (<í>n I ta thu được biểu thức phần bổ chính bậc hai cùa năng lượng: E l ^ X K l H . K K ’ Thay c'1.' = Í Ĩ e ỉÍI I Ĩ ílI vào ta đươc: nk C (0) W 0 ) • E« EL K |H , |O n)f E <2> = y I 11 c(°) _ c(°) (1-75) k*n *--k Để cho phép gần đúng có nghĩa thì số hạng bổ chính phải nhỏ, ta phải có: lc(1)| -K I h .Ịcd,,) 1 Ilk 1H<0’ - E‘k0)« 1 (1.76) 28 5.2. Nhiễu loạn suy biến Xét toán tử Hamilton Ho không nhiễu loạn bị suy biến bậc s H o|® nk) = El01|O nk) , k = l , 2 .... ,s. (1.77) ứng với một giá trị năng lượng E^0), trạng thái cùa hệ không nhiễu loạn được mô tả bời các vectơ trạng thái trực giao: l ^ n Ạ I ^ ) . .....| « 0 Giả sử Hamiltonian nhiều loạn H của hệ có dạng: H = Ho + ẰHi và các vectơ riêng cùa H được xác định từ phương trình Schrodinger. H |T n) = (H o + m ,)|V „ ) = E „ |'ỉ'n) (1.78) Khai triển I'ỉ'n) theo các vectơ riêng cùa toán từ Ho l^ n > = i c k|0>nk> k=l và thay (1.79) vào (1.78) ta được: (1.79) X ck (Ho I ®„k) + X.H, I 4>„k)) = En ị c k I O llk) (1.80) k=I k=l Nhân hai vế của phương trình (1.80) với bra (nk I (* « I ■®„„) +1 « I ^H | I ) = E„ ỉ ; ct (®„, I ®„k) (1.81) k=l k»l Và kí hiệu yếu tố ma trận = (®li; I^.Hi | sẽ tách thành một dãy các mức gần nhau: En = E(n0) + enj , j = 1, 2.....s. Thay một giá trị £ nj vào hệ phương trình (1.82) ta tìm được một giá trị của ck. Thành thử vectơ trạng thái của hệ nhiễu loạn bây giờ có dạng: (1.84) k = l 5.3. Nhiễu loạn phụ thuộc vào thời gian Trong trường hợp nhiễu loạn tác dụng lên hệ lượng tử phụ thuộc vào thời gian. Hamiltonian phụ thuộc vào thời gian của hệ có dạng: H(t) = Ho + Hi(t) (1.85) Năng lượng cùa hệ nhiễu loạn nói chung khòng bảo toàn và hệ không có trạng thái dừng. Tuy nhiên để áp dụng phương pháp nhiễu loạn chúng ta sẽ xác định các vectơ trạng thái của hệ nhiễu loạn theo các trạng thái dừng của hệ khòng nhiễu loạn |n). Trước hết, ta khai triển các vectơ riêng l'i'(t)) cùa Hamiltonian toàn phần theo các trạng thái dừng riêng |O n) cùa toán tử Ho M = Ic„(t)|cDn} ( 1.86) Khi H| = 0 ta có: (1.87) trong đó E„ là năng lượng của hệ ờ trạng thái dìmg |O n) H „ |0 11) = E „ |0 n) 30 Do đó để thuận lợi hơn ta viết lại hệ số c n(t) khi hệ có nhiễu loạn: M o ) = I c n( t) .e x p ^ - ^ j « I > n) ( 1.88) Thay khai triển vào phương trình Schrodinger: '\.h— I 'f (t)) = (Ho + íìi)|'i'(t)) ta được: yidca(t) dt " " ’ ~ ^ l® > l[ E » + H i ( t) ] c n( t) .e ^ |® n) (1.89) Nhân vò hướng (1.89) với bra (O m I và sử dụng tính chất trực chuẩn cùa bracket |n) = Smn, ta được: í 1-9 0 ) Giả sử khi t < 0 hệ ờ trạng thái không nhiễu loạn nào đó |O j) , khi dó: cn(0) = c(nO) =sni (1.91) Bắt đẩu từ thời điểm t = 0 hệ chịu tác dụng của nhiễu loạn nhỏ Hi và ờ thời điểm t > 0 các hệ số cn(t) có thể khai triển dưới dạng chuỗi nhiễu loạn sau: c „ = c („0) + c (n1) +cí12)+... (1.92) Thay (1.92), (1.93) vào (1.90) ta nhận được phương trình đối với cn(t) trong gẩn đúng bậc một ( m * i ). i.(Ein -E j).t i.s < w = exp dt nghiệm của phương trình này có dạng: • t i.(E m- E , ).t' ^ ’ = —y Jdt'e * (Om|Hi (t')|i) do đó trong gần đúng bậc một, ta có: c„,(t) = ômi - J J d t'e " r" ftEi>l (O m |H .(t')|0 ,) với điều kiện: c (t) « 1. (1.93) (1.94) 31 5.4. Phương pháp biến phân Trong cơ học lượng tử, ngoài lí thuyết nhiễu loạn, người ta còn sứ dụng phương pháp gần đúng khác gọi là phương pháp biến phân. Phương pháp biến phân xuất phát từ nhận xét đcm giản rằng năng lượng trung bình cùa một hệ luôn lớn hơn hoặc bằng năng lượng trạng thái cơ bản của hệ lượng từ. Việc tính nàng lượng trạng thái cơ bản dẫn đến việc chọn các hàm thử chứa một sò thõng sô chưa biết nào đó. Sau đó tìm cực tiêu của năng lượng trung bình cho phép ta xác định được thông số, nghĩa là xác định được năng lượng trạng thái cơ bản cúa hệ. Khai triển vectơ trạng thái của hệ lượng tử I 'p ) theo các vectơ riêng I u n) của toán tử Hamiltonian H : |y> = 2 » „ ) (1.95) n Trị trung bình của năng lượng của hệ ờ trạng thái đã cho có dạng:rạng của hệ ( /G \ ( J Ể M . ĩ |a"| : 'E " \ / / v l v ì ^ T O X k l 2(1.96) Gọi E0 là năng lượng trạng thái cơ bản cùa hệ lượng tử, từ (1.96) ta có: Eo Z h f ( h ) > ■ = E 0 (1.97) 2 > . .l n Chọn các vectơ trạng thái Ị1? ) là một hàm cùa thông số chưa biết nào đó X|,sao cho gần trùng với vectơ trạng thái cơ bán cùa hệ. = (1.98) và thực hiện cực tiếu hoá nàng lượng trung bình: c(h\ -± -ỉ- = 0 (1.99) CẰị cho phép xác định các thông số X.;. Bàng cách đó ta sẽ tính được giá trị trung bình . . r- < m „ 1,x,e....)|íì|'í'<>.... x„......) , ... năng lượng E = ——------—----- iy-i-------!— Ỉ— L gần với giá trị Iiãna lượiis cơ bàn ( y của hệ. 32 Phương pháp tính năng lượng ở trạng thái cơ bản của hệ lượng tử nói trên phụ thuộc vào việc chọn hàm thử. Ngoài ra ta cũng có thể tính năng lượng của trạng thái kích thích thứ nhất E, hoặc thứ hai E„... Thực vậy, nếu kí hiệu I '-P(| \ là vectơ trạng thái cơ bản của hệ thì việc tính E, đòi hỏi bài toán biến phân: E, = m in('F1|H |4 /l) (1.100) với điều kiện bổ sung (VF ]|VỈ/1)= 1 v à ( lF1|lF0} = 0 . Tương tự, việc tính năng lượng ờ mức kích thích E2 dẫn đến giải bài toán: E, = m in(4/2|H |T 2) (1.101) với điều kiện bổ sung ( ^ 21^ 2} = 1 và ( ^ 21 ) = ( ^ 21^ 0) = 0 5.5. Phương pháp trường tự hợp Hartee - Fock Toán tử Hamilton của nguyên tử nhiéu electron hệ toạ độ gắn với hạt nhân có dạng: H = Ị H ' + ị r k V u i, k = 1 , 2 ,..., z ( 1 . 1 0 2 ) trong đó H là toán tử Ham ilton của electron thứ i trong trường hạt nhân diện tích e 2 . .... .. . , , ' Ze, v k; = — là toán tử tương tác của hai electron, dấu phay ở tống thứ hai kí hiệu r ki tổng chỉ lấy các giá trị k và i với k * i. Để thuận tiện, ta dùng phương pháp biến phân đê tính năng lượng của trạng thái cơ bản của nguyên tử. Khi đó hàm sóng của nguyên từ được xác định từ đẳng thức. 6J = 5 j\|/*Hv|/dV = 0 (1.103) Với điểu kiện j\ịí*Y|/dV = 1 Ta chọn hàm sóng cùa nguyên tử là tích cùa hàm sóng cùa các electron riêng biệt: V O Ì Â - r2) = (p(r,) ‘kVkii(ri) = 0 (1.111) Nhờ hệ phương trình này, ta tính được thế nãng. Thế năng nhận được như vậy gọi là trường tự hợp Hartee. Như vậy, với phương pháp trường tự hợp, bài toán nhiều electron rút về bài toán một electron. Mỗi electron trong nguyên tử được xét như chuyển động trong trường tự hợp tạo bởi hạt nhân cùng với tất cả các electron còn lại. Fock đã mờ rộng phương pháp của Hartee bằng cách kể đến tính đồng nhất của electron. Theo phương pháp của Fock, hàm thừ được xây dựng nhờ các hàm sóng của electron riêng biệt phụ thuộc cả vào các biến số không gian lẫn các biến số spin. Đặt ^1 là tập V i ( 4 i ) V i ( 4 2 ) .......... v , ( 5 z ) ........V 2 ( 5 z ) \ụ =1 4 Ĩ .( 1.112) v z(4 l)v z(^2)....V z^ z) Fock dã nhận dược phương trình đối với hàm sóng V)/ (q ). Phương trình này có dạng phức tạp hơn nhiều so với phương trình nhận được bằng phương pháp Hartee nhưng kết quả chính xác hơn nhiều. Ví dụ trường hợp nguyên từ Hêli: Tính đối xứng của hàm tọa độ phụ thuộc vào spin toàn phần của hệ. Trong trường hợp nguyên tử Hêli, ta có hai trường hợp: + Trạng thái para ứng với spin toàn phần cùa hệ bằng khòng. + Trạng thái octo ứng với spin tổng cộng bằng một. Giả sử toán từ Hamiltonian có dạng: H = Hi+H°+V |2 (1.113) ~0 Trong đó Hi chi tác dụng lẽn toạ độ cùa hạt thứ ị. 35 — Đầu tiên ta tìm phương trình xác định có trạng thái para của hệ: Trạng thái spin tổng cộng bằng không và hàm sóng tọa độ là đối xứng. Ta chọn hàm thư: y = - ^ k ( l ) < p b(2) + cpa(2) Thay vào biểu thức của J ta có: J= Jvị/* Hvị/.dV, .dV2 = I^ H ÌV d V , + |cp*b.H°. b(2)-íp.,(2)ọb(l)] (1.117) Tương tự, hệ phương trình Fock có dạng: (H + vbb “ E., )cpa - vKltph = 0 + v , a - E h ) chính là momen cấp một còn momen trung tâm cấp hai <(q, - )2> chính là phương sai. Từ (1.119) và (1.120) ta thấy khi biết hàm phàn bô f(q,, q2,..., qn) thì sẽ xác định được các momen. Đối với các hệ lượng tử được mô tả bời các toán tử thống kê P , các momen được định nghĩa như sau: ( 1. 121) ( 1.122) Toán tử P tuân (hco phương trình Liouville lượng tử: r ô ^ = ÍH ,p ] (1.123) ỡt L J tronể đó dấu [...] ờ vế phải là dấu móc Poisson lượng tử. Như vậy, nếu biết toán tử thống kê P thì có thể tính được momen. Mặt khác, giữa các momen có mối quan hệ với nhau, momen cấp cao có thể biểu diễn qua momen cấp thấp hơn. Các hệ thức liên hệ giữa các momen đóng vai trò quan trọng trong việc nghiên cứu các tính chất nhiệt động cùa tinh thể phi tuyến nên phương pháp momen có thể coi là một trong những phương pháp hữu dụng để nghiên cứu bài toán này. 40 PHẦN II NGÔN NGỮ LẬP TRÌNH MATLAB VÀ FORTRAN Chương 2 GIỚI THIỆU VE MATLAB §1. M ỏ ĐẦU Việc sử dụng ngôn ngữ lập trình Matlab mô phỏng một số hệ vật lí là rất hữu dụng, vì đây là ngôn ngữ có khả năng ứng dụng rất lớn và linh hoạt trong quá trình thiết kế. Đé tiếp cận được ngôn ngữ lập trình này trưóc hết chúng ta cần nắm được khả năng ứng dụng, vai trò của ngôn ngữ lập trình trong ngành khoa học kĩ thuật nói chung và ngành Vật lí học nói riêng. 1.1. Matlab - ngôn ngữ của tính toán kĩ thuật MATLAB là một ngôn ngữ bậc cao có môi trường tương tác cho phép tiến hành các nhiệm vụ tính toán có cường độ lớn, nhanh hơn so với các ngôn ngữ lập trình khác như c, C++ và Fortran. MATLAB viết tắt từ "Matrix Laboratory". Ban đầu Matlab dược thiết kế bời ơeve Moler vào những năm 1970, và được sử dụng như một công cụ dạy học. Từ đó đến nay nó đã được phát triển thành một bộ phần mềm thương mại rất thành công. Hiện nay MATLAB R14 là một bộ phần mềm dùng để tính toán trong các ngành khọa học kĩ thuật và toán học ứng dụng. Matlab là một ngôn ngữ lập trình mạnh, có giao diện đồ họa đẹp. Matlab là một thương hiệu đã được thương mại hóa của tập đoàn MathWorks, Massachusetts, USA (hiện là nhà cung cấp hàng đầu thế giới cho các phần mềm tính toán kĩ thuật và thiết kế dựa trên mô hình). 1.2. Khả năng và những ứng dụng của Matlab Một trong những tính năng tuyệt vòi nhất của Matlab nhìn từ góc độ những nhà khoa học tính toán là thư viện dựng sẵn rất phong phú, có các chu trình tính toán và các công cụ hiển thị dồ họa. 41 Matlab cho phép người sử dụng tiến hành rất nhiều các nhiệm vụ thõng thường liên quan tói việc giải quyết các vẫn đề một cách số học. Nó cho phép người sử dụng dành nhiều thời gian hơn cho việc suy nghĩ, khuyến khích làm thí nghiệm. Các tính toán mạnh có thể được thực hiện chỉ với một hoặc hai câu lệnh. D dàng xây dụng riêng những hàm toán học cho những ứng dụng đặc biệt Matlab cung cấp giao diện đổ họa đẹp, các hình từ Matlab có thể chèn vào LATEX và các tài liệu Word. Lập trình theo nghĩa thông thường là nhập vào máy nhũng câu lệnh rõ ràng, theo một thứ tự nhất định sao cho khi máy thực hiện theo đúng thứ lự đó thì sẽ cho kết quả mong muốn. Một khái niệm tương tự như vậy thường thấy trong các ngôn ngữ lập trình khác như: ngôn ngữ c, Pascal... Khi khởi đầu với MatLab ta hãy hiểu theo nghĩa rộng hơn: Lập trình còn có các bước biểu diễn bài toán dưới dạng các hàm và máy tính qua việc thực hiện các hàm này cho ta kết quả. Phương pháp này có mức độ trừu tượng cao hơn so với các câu lệnh chỉ dẫn đơn thuần. Tuy vậy, các ngôn ngữ lập trình biên dịch như c hay Fortran cho phép tính toán nhanh và tốc độ cũng là yêu cẩu quan trọng trong các chương trình tính lớn. Do đó cách kết hợp thông minh là điểu cần thiết phẩn lõi tính toán được viết bằng ngôn ngữ biên dịch, các thao tác nhập xuất, xử lí, hiển thị số liệu được viết bời ngôn ngữ MatLab. §2. CÀI ĐẶT VÀ KHỎI ĐỘNG MATLAB 7.0 2.1. Cài dặt Matlab 7.0 Yêu cầu về cấu hình máy tính: + Bộ vi xử lí Pentium hoặc Pentium Pro. + Windows 95 hoặc NT (WinXP home, XPprofessional). + Bộ điều phối đổ họa 8 bit và card màn hình tôi thiểu 256 màu. + Dung lượng ổ cứng 25Mb cho tới hơn 1Gb (tùy thuộc vào cấu hình đĩa cứng, phân vùng đĩa, số hợp phần của Matlab được cài đặt), và tới 2.1Gb nếu cài đặt Matlab cùng với Simulink. + Bộ nhớ động (RAM) tối thiểu 16Mb (nên có bộ nhớ tối thiểu 128Mb). + Các khuyến nghị khác: Bộ nhớ bổ sung, card đồ họa bo sung, card âm thanh, máy in, MS-Word 7.0 hoặc hơn, trình biên dịch c, Borlean. Microsoft (xây dựng file MEX), trình duyệt internet (đê chạy Matlab Helpdesk online). 42 Quá trình cài đặt Matlab 7.0 cho WindowsXP (bộ gồm 2 đĩa CD): + Đưa đĩa CD vào ổ đọc. Nếu chương trình SETUP không tự động chạy thì nhấn đúp vào biểu tượng setup.exe để bắt đẩu quá trình cài đặt. + Accept (chấp nhận) những thỏa thuận về bản quyền. Sau đó click Next. + Nếu cài theo kiểu mặc định (hay còn gọi là Typical setup - kiểu phổ biến), Matlab trên máy tính sẽ có các hợp phần cơ bản nhất để làm việc theo các hướng dẫn trong tài liệu này. Theo các hướng dẫn trên màn hình. Cho đĩa CD 2 vào khi dược yêu cầu. + Nếu cài đặt theo kiểu tùy chọn cá nhân (Manual setup) thì nhấn vào các hộp thành phần dấu 'v' nếu muốn có tùy chọn đó. Nhấn tiếp nếu không có ý định (có thể thêm vào sau này nếu muốn). + Trên màn hình hiển thị 'C:\MATLAB7' là thư mục mặc định của quá trình cài đặt. Nếu muốn cài đặt vào địa chỉ khác, hoặc đổi tên thư mục, thì ta lựa chọn 'Browse'. + Chi tiết hướng dẫn cài đặt xin xem file ‘install-guide.pdf’ trong đĩa CD1. 2.2. Khỏi động Matlab (Hệ điều hành Windows) Từ hệ điều hành Windows, khởi động Matlab đơn giản bầng cách nháy đúp vào biểu tượng MATLAB trên màn hình, hoặc bằng cách chọn MATLAB từ Menu Start. Quá trình khởi động dưa người dùng đến cửa sổ lệnh, nơi các dòng lệnh được biểu thị bằng dấu Đây là dấu hiệu cho thấy Matlab đang chờ một (câu) lệnh. Khi hoạt dộng trong chế độ máy tính, tất cả các lệnh của Matlab được nhập vào dòng lệnh từ bàn phím. Matlab có thể được sử dụng theo nhiều chế độ và nhiều cách khác nhau: + Như một máy tính cá nhân cao cấp trong chế độ máy tính cá nhân. + Như một ngôn ngữ lập trình bậc cao. + Như một chu trình con gọi lừ chương trình c. §3. QUÀN LÍ KHÔNG GIAN LÀM VIỆC CỦA MATLAB v ể cơ bản, không gian làm việc của Matlab gồm các phần sau: + Cửa sổ trợ giúp (Help window) + Nút Start + Cửa so nhập lệnh (Command window) + Cứa sổ không gian làm việc (Workspace window) + Cửa sổ quá trình lệnh (Command History window - lịch sử) 43 + Cửa sổ biên tập mảng, vectơ, ma trận (Array editor window) + Cửa sổ địa chỉ thư mục hiện thời (Current directory window) - Nút ‘x’ ờ góc trên bên phải mỏi cửa sổ dùng để đóng chúng. Hiến thị lại cửa sổ bằng cách tích vào tên cửa sổ tương ứng trong m enu Desktop. - Nút mũi tên cong bên cạnh nút ‘x ’ dùng dể tách các cửa sổ làm việc trong cửa sổ chính Matlab thành cửa sổ con dộc lập. Ấn nút này một lần nữa sẽ nhập một cửa sổ độc lập về cửa sổ chính của MATLAB. - Nút Start: ở góc dưới bên trái cùa màn hình Matlab, cho phép chạy các ứng dụng mẫu (demos), các công cụ và cửa sổ chưa hiển thị khi khời động Matlab. Bằng cách gõ lệnh 'demo' có thể xuất hiện một tập hợp những file trình diễn giá trị cao. Ví dụ: Thử chạy Start —» Matlab —» Demos để chạy một ứng dụng mẫu trong cửa sổ Demo(s). . MATUD Ht Edl Cwtug OilMp unnoow Htto fItra c *» 1LAB7iw t V (t) D ãỉ X. % m - - ; u cf ý c* SKrtots >]MowloAdd 2J .MX'tMtw f t a 0 <. * * *5 • CurrM t*M»ory yVoilưpoca - 9- iv IT - lr elc i Hình 2.1. Giao diện của Matlah 7.0 Glii chú: S s-D + Các diễn giải và câu (mệnh đề) của Matlab được đánh giá khi gõ vào 'cửa sổ lệnh', và các kết quả tính toán cũng được thể hiện tại đây. Không giống như Fortran và các ngôn ngữ tính toán cẳn biên dịch khác, Matlab là một mõi trường lương tác - đưa ra m ột câu lệnh, thì Matlab cô gắng thực thi câu lệnh đó ngay lập tức trước khi đòi hỏi một lệnh tiếp theo. + Cấu trúc có dạng: » biên = diễn giải hoặc đơn giàn là » diễn giải - Các diễn giải thường được soạn bằng các toán tử, các hàm, tên các biến và được hiển thị trên màn hình sau khi ấn Enter. 44 Các câu lệnh có dạng ‘tên biến = diễn giải’ thì diễn giải đó sẽ được gán cho biến để sử dụng sau này. Khi ‘tên biến’ và dấu ‘=’ được bỏ đi thì kít quả cùa diễn giải sẽ được tự động gán cho biến có tên ‘ans’ (câu trả lời) và hiển thị trên màn hình. Một câu (lệnh) thông thường sẽ kết thúc ở cuối dòng. Tuy nhiên có thể tiếp tục một câu bằng ba dấu chấm ờ cuối dòng. Có thể đặt một vài câu lệnh trên cùng một hàng, ngăn cách bời dấu phẩy Y hoặc chấm phẩy ‘; \ Nếu một câu lệnh kết thúc bằng dấu chấm phẩy ờ cuối câu thì kết quả cùa lệnh đó sẽ không được hiển thị, tuy nhiên yêu cầu tính vẫn được thực hiện (phép tính hay phép gán vẫn dược thực hiện, kết quả có trong workspace). Điều này là thiết yếu trong việc ẩn đi các kết quả trung gian không mong muốn. Có thể xóa trắng toàn bộ cửa sổ lệnh bằng lệnh » clc (clear command window) hoặc vào menu Edit —» Clear Command Window. Khi thực hiện lệnh này, toàn bộ giá trị của các biến hiện có không thay đổi hay mất di. Cửa sổ không giun làm việc (workspace): Các biến và dữ liệu khi nhập vào hoặc sau khi tính toán sẽ được Matlab lưu trong phẩn “không gian làm việc”. Tất cả các biến, ngoại trừ những biến cục bộ thuộc về một M-file, sẽ được thể hiện trong không gian làm việc. - Lệnh 'who' hoặc 'whos' liệt kê các biến hiện có trong không gian làm việc. Để biết giá trị hiện tại của một biến, ta đánh vào tên biến tại dấu nhắc của cửa sổ lệnh và Enter. Để xóa một hàm hoặc biến khòi không gian làm việc, ta sử dụng lệnh 'clear': » clear tên-biến Bản thân lệnh 'clear' sẽ xóa tất cả các biến hiện có (tương đương vói 'clear all') - Cửa sổ biên tập mảng (ma trận nói chung): Khi đã có một mảng, ta có thể chỉnh sửa, biên tập lại bằng Array Editor. Công cụ này làm việc như một bảng tính (spreadsheet) cho m a trận. Cũng có thể biên tập lại ma trận M bằng cách gõ lệnh: » openvar ('C') - Cửa sổ địa chỉ thư mục hiện thời: Thư mục hiện thời là nơi chương trình Matlab sẽ tìm các M-file và các file không gian làm việc (.mat files) mà ta dã Load và Save. a 0 ‘i Ế I T " ta “í ta ì# ^ - ĩtacK ị -(*»•• Ị o*:ĩ M.rf. ir!êtie [1 - 1 4 í.-' 7-. JauHe • \ | ! ............ > Hình 2.2. Cứu sổ không iỊÌan lùm việc 45 §4. MỘT SÔ' LƯU Ý KHI LÀM VIỆC VỚI MATLAB 4.1. Lưu và phục hồi dữ liệu - Để nhớ các biến, Matlab có thể ghi và gọi lại dữ liệu từ file trong máy tính. Mục Save Workplace as... trong bảng chọn File sẽ mờ hộp hội thoại để ghi tất cả các biến hiện tại. - Tương tự, mục Load Workplace trong bảng chọn File sẽ mờ hộp hội thoại để gọi lại tất cả các biến đã ghi lại từ không gian làm việc trước. Ghi chú: Việc Load không làm mất các biến hiện có trong không gian làm việc hiện tại. Khi gọi lại các biến trùng tên với các biến trong không gian làm việc của Matlab, sẽ làm thay đổi giá trị của các biến theo giá trị của các biến gọi ra từ file. Ngoài các bảng chọn, Matlab còn cung cấp hai lệnh Save, Load nó thực hiện một cách mềm dẻo hơn. Lệnh save cho phép ghi một hoặc nhiều hơn một biến tùy theo sự lựa chọn. Ví dụ: » save - lưu tất cả các biến trong Matlab theo kiểu nhị phân trong file matlab.mat » save dulieu — lưu tất cả các biến trong Matlab Iheo kiểu nhị phân trong file dulieu.mat » save dulieu A B c D - ASCII lưu các biến A, B, c , D theo dạng mã ASCII trong file dulieu.mat. - Lưu lại một phiên làm việc (session) Khi làm việc vói Matlab, cần lưu tất cả các thông số đầu vào và đầu ra của phiên làm việc để thuận lợi cho việc in ấn sau này. Lệnh ’diary’ sẽ lưu tất cả những thông số đầu vào và đẩu ra ờ giữa hai lệnh 'diary' và 'diary off. Ví dụ: » diary('diary-file-nam e') » .........(các câu lệnh của bạn ờ đây) » diary off 4.2. sửdụng Help (Trợ giúp) Trợ giúp thõng tin về các lệnh cùa Matlab có thê được tìm thấy theo nhiều cách: + Từ dòng lệnh bằng cách gõ lệnh help chủ đề'. 46 + Từ cửa sổ màn hình vào Menu Help. + Từ helpdesk của Matlab lưu trữ trên đĩa hoặc CD-rom. + Từ dòng lệnh gõ lệnh 'help' và Enter. Kết quả: Matlab cho một bản tóm tắt về hệ thống trợ giúp HELP topics: (Các chủ đề trợ giúp) matlab/gẽneral - Các lệnh với mục đích tổng quát, matlab/ops - (operators) Các toán tử và các kí tự đặc biệt... matlab/lang - (language) Ngôn ngữ lập trình... matlab/elmat - (elementary) Ma trận căn bản... matlab/elfun - (elementary functions) Các hàm toán căn bản. matlab/specfun - (specialized functions) Các hàm toán đặc biệt. - Thõng thường cửa sổ text không đủ lớn để chứa tất cả thông tin từ một lệnh Matlab. Vì vậy có thể sử dụng chức năng 'more on' để xem thông tin theo từng trang màn hình, sau đó duyệt từng trang một bằng cách nhấn phím bất kì. - Đánh 'more off vào cửa sổ lệnh sẽ đưa Matlab trở về cửa sổ thông thường, không duyệt từng trang. Khi không nhớ chính xác tên của một lệnh Matlab, có thể sử dụng lệnh 'lookíor' (tìm kiếm) như một sự trợ giúp. ■> Help B E ® 1 Ffe Ed* View Go Favw*« Desktop Window H<* Hfelp Mđvigcípr Xf 4 C â N Coríerts irvđey í Search Oemos __________ Titt*: raatog (MATLAB Functons) V Seereỉì for. logafOvn V Go ......... ........ ^ A| Tile I Section MATLAB Function Reference ® f f l 1 Example: Elementwise Logarithm Galois Field Computations A 1 r e a l l o g betaln MA|TIAB Functions Linear Algebta Functions -- Categoncal Lis j Natural logarithm for nonnegatwe real arrays Elementary Math Functions -- Categoncal Lis’ 1 Specialized Malh Functions ” Categorical Lis S y n ta x funm MATLAB Functions Y = r e a llo g ( X ) gallery MATLAB Functions gamma, gammainc, gammaln MATLAB Functions D e s c rip tio n gamma, gammainc, gammaln MATLAB Function? gamrna. gammainc, gammaln MATLAB Functions Y “ re a llo g C S ) returns the natural log MATLAB Functions logarithm of each etenent in array X Array X logio MATLAB Functions must contain only rionriegaitve real numbers log2 MATLAB Functions The 5Ĩ28 of Y is the same as the size of X login MATLAB Functions V 1 < > j E x a m p le s £S3ICtLSjMlj2£lJ29Mi5S.to W6b If r M = m agic (*i) V looaiithm < * 43 pages contain the vvc-rd logarithm Hình 2.3. Giao diện cửu sổ Help của Matlab 7.0 47 4.3. History Trong giao diện của Matlab, cửa sổ 'Command History' nằm ờ góc phần tu bèn dưới, phía trái. Trong cửa sổ này, các lệnh đã sử dụng trong các lần khởi động Matlab gần đây đều được lưu lại. Mỗi lần khời động Matlab, toàn bộ các lệnh sừ dụng trong lẩn đó sẽ được lưu lại dưới dạng rr.ột nhóm có thể đóng mớ bằng nút biếu tượng '+' (expand) hoặc (collapse) ở đầu dòng (dòng ghi m ốc thòi gian giữa hai dấu chú thích '%'): Ví dụ: (+) % - 8/14/07 3:21 PM -% Để gọi lại lệnh từ cửa sổ 'Command History’, bạn tìm đến lệnh đó bầng các thanh cuốn, rồi nháy đúp vào tên lệnh. Đê gọi lại các lệnh bạn đã sừ dụng từ đấu nhắc của cửa sổ lệnh, M atlab dùng các phím mũi tên (<-T—>-L ) trên bàn phím. 4X Chương 3 NGÔN NGỮ LẬP TRÌNH TÍNH TOÁN TRÊN MATLAB §1. CÁC PHÉP TOÁN Cơ BẢN Trước hết chúng ta cần nắm bắt được cách sử dụng các phép toán với các biến, các phép toán có thể từ những phép toán sô' học thông thường (Phép cộng, trừ...) đến các phép toán phức tạp như: Giải phương trình vi phân... còn các biến có thể bằng sô' hoặc chữ. 1.1. Các phép toán số học thông dụng: Cộng (+), trừ (-), nhàn (*), chia (/) Ví dụ 1: Để tính A=2+3/4*5 B=2A5 - 3*A A+B Bạn gõ các lệnh vào sau dấu nhắc lệnh: » - » A = 2 + 3/4*5 A = 5.7500 » B = 2A5 - 3*A B = 14.7500 » A + B ans = 20.5000 Chú ý: - Khi không có toán từ gán, Matlab trả kết quả của phép tính gần nhất vào biến 'ans= - Thứ tự ưu tiên tính toán: Các đại lượng trong ngoặc đơn, lũy thừa, (* /) làm việc từ trái qua phải, (+ - ) làm việc từ trái qua phải. Lệnh gán: - Gán trực tiếp các giá trị cho các biến - Ngoài ra: Có thể dùng biến trả dế gán Ví dụ 2: » 2A3-3 ans = 49 5 » ans*6 ans = 30 Quy tắc đặt tên biến: - Tên biến hợp lệ được cấu tạo bới các chữ và số, bắt đầu bằng chữ. - Chiéu dài tên biến: Độ dài tùy ý Chú ý: Tên biến không phụ thuộc vào chữ hoa - chữ thường Tránh đặt tên biến trùng với tên các hàm chuẩn, hoặc các từ khóa cùa Matlab như: 'break' 'case' 'catch' 'continue' 'else' 'elseif 'end' 'for' function' 'global' 'if 'otherwise' 'persistent' 'return' 'switch' try' 'while'... 1.2. Các hàm số cơ bản Lũy thừa: A, giai thừa: x! Logarit: ln(x), log[a](b), exp(x) Các hàm lượng giác: sin(x), cos(x), tan(x), cot(x),... Một số hàm khác: abs(x) - Ixl, sqrt(x) - căn bậc 2 cùa X 1.3. Các phép toán quan hệ 1.3.1. Các toán tử so sánh A < B % A nhỏ hem B A <= B % A nhỏ hem hoặc bằng B A > B % A lớn hơn B A >= B % A lớn hơn hoặc bằng B A == B % A bằng B A ~= B % A khòng bằng B - Các toán tử quan hệ Ihực hiện sự so sánh từng phần tử với phần từ giữa hai máng. Nó cho kết quà là một máng logic có cùng kích cỡ với các phần tứ cùa mảng là đúng (l) nếu quan hệ đó là đúng, và phán từ cúa mảng là sai (0 ) nếu không đúng. 1.3.2. Các phép toán logic - Matlab biểu diễn đúng và sai bời các số nguyên tố l và 0: đúng = l. sai = 0 (lrue=l, false=0). Ví dụ 1: Trong quá trình tính toán, biến X (x là một đại lượng vò hướng) nhận một giá trị bất kì. chúng ta có thế tiến hành các phép kiêm tra logic ch o nó: 50 X = = 2 xem X có bằng 2 không? X ~ = 2 xem X có khác 2 không? X > 2 xem X có lớn hơn 2 không? X < 2 xem X có nhỏ hơn 2 không? X >= 2 xem X có lớn hơn hoặc bằng 2 không? X <= 2 xem X có nhỏ hơn hoặc bằng 2 không? 1.3.3. Các toán tử logic chinh trong Matlab và ý nghĩa của chúng & và (and) hoặc (or) không (not) 1.3.4. Các kiểu dữliệu - Định dạng kết quả 1.3.4.1. Các kiểu dữ liệu Matlab sử dụng 15 kiểu (loại) dữ liệu chính. Mỗi một kiểu dữ liệu đều ở dạng của một ma trận hoặc mảng. Các mảng hoặc ma trận có kích cỡ tối thiếu là 0 - nhân - 0 và có thể phát triển tới mảng n - chiểu với kích cỡ tùy ý. a. Các kiểu dữ liệu số và số phức - Integer: ví dụ như -5 hay 9888. - Double precision reals: Trong Matlab, tất cả các số thực được lưu với độ chính xác double, không giống các ngôn ngữ lập trình khác như c hay Fortran khi chỉ có một loại riêng biệt float hay real cho các số thực với độ chính xác single. - Số phức: Được nhập vào dưới dạng 3+2*i hoặc 3+2*sqrt(-l). - Chuỗi: Là một mảng tập hợp của các kí tự, được nhập vào dưới dạng 'abc' hoặc 'vi du day la mot chuoi'. b. Các kí tự, chuỗi và văn bản Khả năng xử lí văn bản trong tính toán rất hữu ích cho việc nhập (xuất) kết quà từ (tới) màn hình vào file lưu Irên đĩa. Đế có thể quản lí văn bản, một loại dữ liệu 'character' được đưa vào Matlab. Một mảng đơn giản của văn bản là một chuồi (vectơ) hay một mảng các kí tự. Các chuỗi có thế được cộng với nhau bằng cách sử dạng các toán tử thao tác trong mảng. Chú ý: Số các kí tự ờ hai dòng phái bằng nhau, nếu khỏng việc thực thi câu lệnh sẽ dẫn tới lỗi: Dấu ba chấm thê hiện câu lệnh còn tiếp tục ớ dòng sau. 51 c. Chuyển đổi giữa chuỗi và số Đôi khi chúng ta cần chuyển một chuỗi thành một số tương ứng hoặc ngược lại. Các công việc chuyển đổi này được Ihực hiện bời lệnh: - 'str2num': Chuyển một chuỗi thành số tương ứng - 'num2str': Chuyển một số thực thành chuỗi tương ứng - 'int2str': Chuyển một chuỗi thành số tương ứng d. Các hằng số Matlab định nghĩa sẩn nhiều hàm số hữu ích, bao gồm: pi = 3.141592654... i và j đều là phần ảo của số phức, = sqrt(-1) inf, 'infinity' là 'vô cùng' NaN, 'not-a-number' là 'không phải là một số' ans luôn được gán cho kết quả cùa lệnh tính trước dó e. Các hàm dựng sẵn Cũng như ngôn ngữ bậc cao khác, Matlab thực thi các 'function' (hàm) nhiều hơn 'procedure' (chương trình con). Các hàm này bao gồm căn bậc hai (sqrt), lũy thừa (exp), logarit (log, log 10, log2), giá trị tuyệt đối (abs) và các hàm lượng giác (sin, cos, tan, cotan,..) Ví dụ 1: » sin(45) ans = 0.8509 Chú ý rằng tất cả các tính toán cùa Matlab đểu có lỗi làm tròn. Danh mục các hàm dựng sẵn phổ biên: + Các hàm lượng giác: + sin - hàm sin, sind - sin cùa góc tính theo độ, sinh - sin hypecbolic, asin - arcsin hay hàm nghịch đáo cùa hàm sin, asind - hàm nghịch đào của hàm sin, kết quả theo độ, asinh - hàm nghịch đáo cùa hàm sin hyperbolic. + COS - hàm cos, cosd - cos của góc tính theo độ, cosh - cos hyperbolic, acos - hàm nghịch đảo của hàm COS, acosd - hàm nghịch đảo của hàm COS, kết quả theo độ, acosh - hàm nghịch đảo của hàm COS hyperbolic, + tan - hàm tang, tand - tang của góc tính theo độ, tanh - tang hyperbolic atan - hàm nghịch đảo của hàm tang, atand - hàm nghịch đảo của hàm tang, kết quả theo độ, atan2 - hàm nghịch đảo của hàm tang 4 góc phần tư, atanh - hàm nghịch đảo của hàm tang hyperbolic, + cot - hàm côtang, cotd - côtang của góc tính theo độ, coth - côtang hyperbolic, acot - hàm nghịch đảo của hàm côtang, acotd - hàm nghịch đảo của hàm côtang, kết quả theo độ, acoth - hàm nghịch đảo của hàm côtang hyperbolic. + Các hàm sơ cấp: exp - hàm mũ, expm l - tính chính xác e x p (x )-l, log - logaril cơ số tự nhiên, loglp - tính chính xác log(l+x), log 10 - logarit cơ sô' 10, reallog - loga cơ sô' tự nhiên của sô' thực, realsqrt - cân bậc hai của một số >= 0, sqrt - căn bậc hai, nthroot - nghiệm thực bậc n của các số thực. + Các hàm liên quan đến số phức: abs - giá trị tuyệt đối, angle - góc pha, 53 + com plex - xây dựng đữ liệu vể số phức từ các phần thực và ào, conj - liên hợp của phức, imag - phẩn ảo của phức, real - phẩn thực cùa phức, isreal - hàm logic, trả về giá trị 'true' với mảng số thực, cplxpair - sắp xếp các số vể các cặp liên hợp phức. + Các hàm làm tròn và phần dư: fix - làm tròn về phía 0, floor - làm tròn về phía âm vô cùng, ceil - làm tròn vé phía dương vô cùng, round - làm tròn vể phía số nguyên gần nhất, mod - mô đun (lấy phần dư của phép chia), rem - lấy phần dư của phép chia (tương tự m od), sign - hàm lấy dấu cùa một biến, trả về +1, 0, - 1 (+, 0, - ) . f. Kết hợp nhiều lệnh trên một dòng Dấu phẩy (,) và dấu chấm phẩy (;) là những kí tự có ý nghĩa đặc biệt trong Matlab như: + Toán tử phẩy (,) được dùng dế nhóm nhiều lệnh trên một dòng. Ví dụ 1: » x = 3.5, y= -5.0, XA 3 - y Khi khôna muốn theo dõi kết quả các tính toán trung gian hay ẩn đi một càu lệnh, ta dùng dấu chấm phẩy (;). Ví dụ 2: » x = 3 . 5 ; y = - 5 . 0 ; X A 3 - y ans = 47.8750 Trong ví dụ trẽn kết quá cùa hai lệnh gán đầu tiên đã được ẩn đi. 1.3.4.2. Định dạng kết quả - Sử dụng lệnh 'format' cùng các định diing. Lệnh này chi làm thay dổi kết quá được hiên thị trên màn hình, khõns làm thay đổi độ chính xác cua sô hoặc phép tính. - Đê thực hiện lệnh, từ díu nhắc cửa sổ lệnh 20 một trong các lệnh sau: format short : chill phẩv Ihập phim cố định, 5 chữ số. 54 format long : dấu phẩy cố định, 15 chữ số, format short e : kí hiệu khoa học, 5 chữ số, format long e : kí hiệu khoa học, 15 chữ số, format short g : dấu phẩy cô' định hoặc di động, 5 chữ số, format long g : dấu phẩy cố định hoặc di động, 15 chữ số, format hex : format dạng Hexa (hệ 16), format V : dương (+), âm (-) và kí tự trắng (blank) ứng với 0, format bank : Dollars và cents, format ra t: tỉ lệ xấp xỉ integer. format short là dạng mặc định, khi được gọi lên, một dạng format sẽ có hiệu lực tói khi nó được thay đổi. 1.4. Giới hạn Các lệnh giới hạn gổm có: - limit(f): giới hạn của hàm f, - lim it(f,x,a) hoặc limit(f,a): giới hạn của hàm f khi X dần tới a, - lim it(f,x,a,’left’): giới hạn trái của hàm f khi X dần tới a, - lim it(f,x,a,’right’): giới hạn phải của hàm f khi X dần tới a. 1.5. Giải phương trình Sừ dụng lệnh solve(f) Ví d ụ 1: + » so ]v e (lf(x) = 0') để giải phương trình f(x) = 0. » so lv e ('2x + 1') » a n s = -1/2 + solve(‘ f(x) = g(x) ’) để giải phương trình f(x) = g(x) » s = solve('cos(2*x) + sin(x) = 1') s = t 0] [ pi] [ l/6*pi] [ 5/6*pi] + solve('f(x)Yg(x)Yh(x)’,...): giải hệ nhiều phương trình. 55 1.6. Tính tích phản + Để tính tích phân ta dùng hàm quad (tính tích phàn theo phương pháp Simpson) và hàm quadl (tính tích phân theo phương pháp Lobatto). Ví dụ 1: f = inỉine('lJ((x-0.3)A2+0.0l)+l./((x-0.9).*2+0.04)-6 '); - q = quuJ(ffi,l) 0) y = prod(l:x-l); else y = gamma(x); end 1.8.2. Câu lệnh switch...case...case... case...otherwise...end Bộ câu lệnh if / elseif chỉ hữu ích trong trường hợp chỉ có một vài lựa chọn. Còn khi có một số lượng lớn các lựa chọn khả đĩ người ta dùng switch để thay thế. Ví dụ: switch donvi case ‘Chieudai’ disp('met’) case ‘The tich’ 57 disp(‘lit’) case ‘Thoi gian’ disp(‘giay’) otherwise disp(‘Toi chiu thua’) end Trường hợp phù hợp với case thì các lệnh sẽ được thực thi. Chú ỷ: Có thể sử dụng otherwise hoặc không. Trong trường hợp có sử dụng thì Matlab thực thi các lệnh sau otherwise, hoặc không nếu không có trường hợp nào phù hợp với case. 1.8.3. Vòng lặp (for và while) 1.8.3.1. Vòng lặp for...end Vòng lặp for., end được sử dụng khi lặp một đoạn mã lệnh với một số lần tùy ý. Ví dụ 1: Vẽ đồ thị hàm số sin(n.x); với -1 < X < 1 và n = 1, 2, . . .8 . Đê thực hiện, ta phải sử dụng 8 lệnh vẽ riêng rẽ. Nhưng sẽ dễ dàng hơn nhiều nếu ta sử dụng một vòng lặp for...end. » X = —1:0.5:1; % bước nhảy 0.5 » for n = 1:8 subplot(4,2,n), plot(x,sin(n*pi*x)) end // \ / \ \ / . \ /J -0 5 ũ 0 5 - - 0 5 0 0.5 1 / \ / \ / \ i \ ị \ / ' V \ \0/\ /\ /\ /\ \ / \ l \ / \ , -0 5 0 Q5 - -0.5 0 0.5 i\ ỉ\ í\ (\ A -/ \j \ / \/ \/ 0'\ A A A !\ A V V \/ ự \ỉ ■ -G 5 0 0 5 - ■0.5 0 0.5 0 ' A /» J ì'1h í'i *' V V V ./ V V V Y " V é ’ i 05 ' 1 -ì ' -0 5 ' Ố Ò T I 58 Chú ý: + Tất cả các lệnh giữa hai dòng bắt đầu bằng ‘for’ và kết thúc bằng ‘end’ đều được lặp đi lặp lại với n = 1 lần thứ nhất, n = 2 lần thứ 2... cho tới khi n - 8. + Lệnh subplot tạo ra một mảng 4x2, cửa sổ đồ thị con trong một đồ thị chính, ở lần lặp thứ n, một hình sẽ được vẽ lên cửa sổ đồ thị con thứ n. 1.8.3.2. Vòng lặp while...end Đẽ thực hiện việc lặp đi lặp lại một đoạn mã lệnh của Matlab cho tới khi một điều kiện (logic) nào đó dược thỏa mãn mà không biết số lẩn lặp cụ thể, khi đó ta có thể sử dụng vòng lặp while...end Ví dụ 1: Tim giá trị lớn nhất cùa n sao cho tổng dưới đây nhỏ hơn 100? » s = 1; n = 1; » while s+ (n+l)A2 < 100 n = n+1; s = s + nA2; end » [n, S] ans = 691 1.9. Vectơ và Ma trận 1.9.1. Vectơ 1.9.1.1. Giới thiệu Đại số tuyến tính là trái tim cùa Matlab. Vì vậy hơn bất kì ngôn ngữ nào khác, Matlab khuyến khích ta tận dụng mọi khả năng của mảng, vectơ và ma trận. * Một vài thuật ngữ: + Mảng là một tập hợp các số, được gọi là các ‘phần tử’ hay các ‘đầu số’, được biết đến với một hoặc nhiều chỉ sô' chạy suốt các tập hợp chỉ số. + Số chiéu của một mảng là sô các chỉ số cần thiết đê định nghĩa một phần từ trong màng. Chẳng hạn mảng 2 chiểu sẽ cần 2 chỉ số i và j đê đặc trưng cho một phần tử của mảng. + Kích thước của mảng là một danh sách các kích thước cùa các tập hợp chi số. V íd ụ l: » r = [ l 2 3; —1 -2 -7 | » size(r) nghĩa là kích Ihước của màng r sẽ là 2x3 (2 hàng, 3 cột). + Ma trận là một máng hai chiều (kích thước m X n với các quy luật đặc biệt cho phép cộng, nhân và các tính toán khác). 59 + Vectơ là một ma trận một chiều có chỉ sô' bằng 1. + Các phần tử đơn lẻ trong ma trận có thể được tiếp cân hoặc sừa đổi bang cách sử dụng chỉ số phần tử (subscripting). Trong Matlab, phần tử thứ i của vectơ V được biểu diễn bằng kí hiệu V(i), chỉ số được viết trong ngoặc dơn. V í dụ 2: » V = [10 20 30] V = 10 20 30 » V(2) ans = 20 1.9.1.2. Vectơ hàng + Vectơ hàng là chuỗi các số được phân cách bời dấu phẩy hoặc khoảng trông. Sô lượng các đẩu sô được gọi là ‘chiều dài’ của vectơ, mỗi đầu sô thường dược nhắc đến như ‘phần tử’, hoặc ‘hợp phần’ của vectơ. + Cú pháp: Để nhập 1 vectơ là một chuỗi các giá trị được bao trong cặp ngoặc vuông [ ]. V íd ụ l: » V = [ 1 3 sqrt(5)] + ngoài ra, ta có thể sử dụng các dấu phẩy (,) » V = [1, 3, sqrt(5)] 1.9.1.3. Vectơ cột Vectơ cột có cấu tạo tương tự như vectơ hàng. Khi định nghĩa vectơ cột, các phần tử được phân cách nhau bởi kí tự hoặc bời ‘newlines’. V í dụ: » c = [ 1; 3; sqrt(5)] 1.9.1.4. Toán tử hai chấm (: ) Toán tử này dùng đê tạo ra vectơ hàng một cách nhanh chóng: » X = 1 : 4 Tổng quát: a : b : c sẽ tạo ra một vectơ vói các phần tử bắt đầu từ giá trị cùa a, tăng dần với bước tăng bằng giá trị của b, cho tới khi đạt tới giá trị của c. 60 » 0 .3 2 :0 .1 :0 .6 ans = 0.3200 0.4200 0.5200 Toán tử còn được dùng để trích xuất một phần của vectơ. Giả thiết chúng ta có vectơ » r = [1:2:6,-l:-2:-7] r = 1 3 5 -1 -3 -5-7 Để trích ra các phần tử từ thứ 3 đến thứ 6 ta có thể dùng lệnh: » r(3:6) ans = 5-1 -3-5 1.9.1.5. Làm việc với vecto và ma trận u. Các phép toán số học + Cộng, trừ ứng với các vectơ có cùng chiều dài. V í dụ 1: » v l =[1 2 3]; » v2 = [4 5 6]; » vl+v2 ans = 5 7 9 + Một vectơ cũng có thể nhân được với một đại lượng vô hưống (một số), bằng cách Matlab tiến hành với từng phần tử. Đê thực hiện các phép tính toán: Nhân, chia và lũy thừa, Matlab đưa ra các toán tử *./’ và \ A\ Ví dụ 2: » vl.*v2 ans = 4 10 18 + Toán tứ lũy thừa có thê được sử dụng Iheo hai cách, lũy thừa sô vô hướng hoặc lũy thừa vectơ: » v2.A2 ans = 61 16 25 36 Trong đó: + Dấu (.*) có ý nghĩa của một phép nhân ma trận hay vectơ. + Tất cả các hàm số học dựng sẵn cùa Matlab được thiết kế để hoạt động vói các vectơ và ma trận, vì vậy chúng ta có thể xây dựng các diễn giải đại sô' hoạt động với từng phần từ của vectơ. Ví dụ 3: » x = [l 2 3];y = [4 5 6]; » s = 2*sqrt(x) + x./y - x.A3.*cos(pi*y) s = 1.2500 11.2284-23.0359 Chú ý: Các phép tính cùa các đại lượng vô hướng trên các vectơ khác với cách làm việc phẩn tử với phần tử. 2*sqrt(x) là nhân số vô hướng với vectơ, trong khi x/y thì ta cẩn phải sử dụng (x./y). + Các phép cộng và trừ phần tử với phần từ không cần thiết phải sử dụng (.+) và h. Ghép các vectơ Đê tạo ra một vectơ từ những vectơ có trước thì kích thước của chúng phải tương thích với nhau. Ví dụ 1: » w = [l 2 3 ],z = [8 9] » cd = [2*z,-w] » so rt(cd ) w = 12 3 z = 8 9 cd = 16 18 -1 - 2 - 3 ans = - 3 - 2 - 1 16 18 62 Trong đó: Câu lệnh cuối cùng (sort) dùng để sắp xếp các phần tử của vectơ theo chiều tăng dần. Ngoài ra, ta cũng có thể sử dụng các lệnh cat, vertcat, horzcat để ghép nối các vectơ. c. Cúc lệnh clio biết thông tin về ma trận (vectơ) size - kích thước theo mỗi chiều length - kích thước của chiều dài nhất (đặc biệt là cho vectơ) ndims - sô' chiều find - các chỉ số cùa các phần tử khác 0 d. Chuyên vị ma trận Để chuyển đổi một vectơ hàng thành một vectơ cột (và ngược lại) bằng một quá trình gọi là ‘chuyển vị’ kí hiệu bằng: (’). Ví dụ 1: » w, w’, c, c' w = 1 -2 3 ans = 1 -2 3 c = 1.0000 3.0000 2.2361 ans = 1 . 0 0 0 0 3 . 0 0 0 0 2 . 2 3 6 1 e. Xiỉrlí dữ liệu với các hàm dựng sun cho vectơ - sort sắp xếp dữ liệu, - max tìm giá trị lớn nhất, - min tìm giá trị nhỏ nhất, - sum tính tổng, - mean tìm giá trị trung bình, - std tìm độ lệch quân phương. 63 1.9.2. Ma trận 1.9.2.1. Cú pháp + Giống với vectơ, các khoảng trông (hoặc dấu phẩy) phân cách các phần tủ trong một hàng và các dấu chấm phẩy là kí hiệu cho biết bắt đầu một hàng mới sau đó. V íd ụ l: » A = [2-1 00 ; 1 I 2 3 ;-1 4 0 5 ] Matlab sẽ đưa ra kết quả A = 2-100 1123 - 1 4 0 5 + Các phần tử đơn ỉé của một ma trận có thể được chỉnh sửa như với các vectơ. V í dụ 2: Lệnh A(3,2) = 0 sẽ thay thế giá trị phần tử cột 2, hàng 3 của ma trận A thành 0. Có nhiều lệnh để khởi tạo một số dạng ma trận đặc biệt zeros(n,m) - tạo ma trận với tất cả các phần tử = 0, ones(n,m) - tạo ma trận với tất cả các phần tử = 1, eye(n) - tạo ma trận đcm vị kích thước n X n. + Để khởi tạo một ma trận vuông đặc biệt, ta có thể sử dụng dạng ngắn zeros(n), lệnh này ngầm định số hàng và số cột của ma trận là bằng nhau. 1.9.2.2. Một số ma trận đặc biệt eye - ma trận đơn vị, zeros - ma trận với tất cả các phần tử = 0, ones - ma trận với tất cà các phần tử = 1, diag - ma trận đường chéo (hoặc chích xuất một đường chéo), toeplitz - ma trận với mỗi đường chéo bằng một hằna số, trĩu - ma trận tam giác trên, tri] - ma trận tam giác dưới, rand - ma trận với các phẩn tử ngẫu nhiên (từ -1 đến 1), linspace - ma trận với các phần tử cách đều nhau, cat - móc nối các ma trận với nhau theo một chiểu đã định, repmat - xây dựng ma trận mới bằng cách lặp một veclơ theo một chiều. 64 1.9.2.3. Các phép toán vói từng phần tử trong ma trận Cộng: x = A + B Trừ: x = A - B Nhân: x = A * B X.*A nhân các phần tử tương ứng với nhau Chia: X = A/B lúc đó X*B = A X = A\B lúc đó A * x = B X = A./B chia các phẩn tử tương ứng với nhau Luỹ thừa: X = AA2 và X = A.A2 Nghịch đảo: X = inv(A) Định thức: d = det(A) cond(A) điều kiện norm(A) chuẩn hoá ma trận eig(A) giá trị riêng, vectơ riêng. §2. CÁC PHƯƠNG PHÁP GẦN ĐÚNG 2.1. Phương pháp chia đôi Để sử dụng phương pháp chia đôi ta xây dựng hàm bisectionO function [x,err,xx] = bisection(f, a, b, tolx, maxiter) %bisection.m de giai pt f(x) = 0 bang phuong phap chia doi cung %vao: f = ham can tim nghiem % a/b = bien cua doan can tim nghiem % tolx = sai so mong muon % maxiter lan lap max %ra: X = nghiem % err = sai so % XX = cac gia tri trung gian tol = eps; fa = feval(f, a); fb = feval(f, b); if fa*fb > 0 error(‘Nghiem khong o trong doan nay’); end for k = 1: maxiter xx(k) = (a + b)/2; fx = feval(f, xx(k)); eư = (b - a)/2; if abs(fx) < tol I abs(err) < tolx break; elseif fx*fa > 0 a = xx(k); fa = fx; else b = xx(k); end end X = xx(k); if k = = maxiter fprintf(‘Khong hoi tu sau %d lan lapSn’, maxiter), else fprintf(‘Hoi tu sau %d lan lap\n’,k), end Ví dụ 1: Để tìm nghiệm của hàm f(x) = tan(pi - x) - X clear all, clc f = inline(‘tan(pi - x) - x ’); [x, ss, xx] = bisection(f, 1.6, 3, le -4 , 50). 2.2. Phương pháp lặp Để sử dụng phương pháp lặp ta dùng hàm simpiterO function [x, err, XX] = simpiter(g, xO, tolx, maxiter) % giai pt X = g(x) tu xO bang cach lap %vao: g, xO = ham va gia tri dau % tolx = sai so mong muon % maxiter = so lan lap max %ra: X = nghiem 66 % err = sai so lx(k) — x(k - 1)1 % XX = cac gia tri trung gian if nargin < 4 maxiter = 100; end if nargin < 3 tolx = le-6; end xx(l) = x0; for k = 2: max iter xx(k) = feval(g, xx(k - 1)); err = abs(xx(k) - xx(k - 1)); if err < tolx break; end end X = xx(k); if k == maxiter fprintf('Khong hoi tu sau %d lan lapW, maxiter) else fprintf^Hoi tu sau %d lan lapNn^k) end Ví dụ 1: Để tính lại ví dụ trên ta dùng hàm simpiterO clear all, clc f=inline('-0.5*((x- l).A2 -3 )'); [x, ss, XX] = sim piter(f, 0.5,.00001,200). 2.3. Đa thức nội suy lagrange Dùng hàm lagrange() để thực hiện việc nội suy hàm số: function [/, L] = lagrange(x, y) %Dua vao: X = [xO xl ... xn], y = [yO y 1 ... yn] %ket qua: / = He so cua da thuc Lagrange bac n % L = Da thuc Lagrange 67 n = length(x) -» 1; %bac cua da thuc/ 1 = 0; form = l: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 (m ,:) = p; %da thuc Lagrange / = 1 + y(m)*p; end Ví dụ 1: Cho hàm dưới dạng bảng: X -2 -1 1 2 Y -6 0 0 6 và tìm y(2.5) ta dùng hàm largrange(): clear all, clc X = [—2 - 1 1 2]; y = [-6 0 0 6]; I = lagrange(x, y); yx = polyval(l, 2.5) 2.4. Đa thức Newton Dùng hàm newton() để nội suy: function [n,DD] = newton(x,y) %Dua vao : X = [xO xl ... xNỊ % y = [yO yl ... yN] %Lay ra: n = he so cua da thuc newton bac N N = length(x) - 1; DD = zeros(N + I, N + 1); DD( 1 :N + 1, 1) = y/ ; for k = 2:N + 1 68 for m = 1: N + 2 - k DD(m,k) = (DD(m + 1, k - 1) - DD(m, k - l))/(x(m + k - 1) - x(m)); end end a = D D (1,:); n = a(N+l); fork = N :-l:l n = [n a(k)] - [0 n*x(k)]; end Ví dụ 1: Cho hàm dưới dạng bảng: X -2 -1 1 2 4 Y 0 0 6 60 Ta dùng hàm newton() để 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) 2.5. Cách tìm cực trị của hàm Phương pháp newton: 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 để cho gradient g(x) của hàm f(x) bằng zero. Nghiệm của g(x) = 0 tìm được bằng cách dùng phương pháp newton cho hệ phương trình phi tuyến. Hàm newtons(x) dùng để tìm nghiệm của phương trình g(x) = 0 là: function [x, fx, xx] = newtons(f, xO, tolx, maxiter) h = le-4; tolfun = eps; EPS = le-6; fx = feval(f, xO); nf = length(fx); nx = length(xO); if nf ~= nx 69 error^Kich thuoc cua g va xO 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]A-l* fx xx(k + 1,:) = xx(k,:) + dx/; fx_= feval(f, xx(k + 1,:)); fxn = norm(fx); if fxn < tolfun I norm(dx) < tolx break; end end X = xx(k + 1,:); if k == maxiter fprintf(;Ket qua tot nhat sau %dlan lap\n', maxiter) end function g = jacob(f, X, h) %Jacobian cua f(x) if nargin < 3 h = le—4; end h2 = 2*h; n = length(x); X = x(:)/; I = eye(n); for 11= 1 :n g(:, n) = (feval(f, X + I(n, :)*h) - feval(f, X - I(n,:)*h))'/h2; end 70 V íd ụ l : Để tìm cực tiểu của hàm bằng phương pháp newton ta có: clear all, clc f = inline('x(l).A2 - x(l)*x(2) - 4*x(l) + x(2).A2 - x(2)'); g = inline(/[(2*x(l) - x(2) -4) (2*x(2) - x (l) - 1)]'); xO = [0.1 0.1]; tolx = le-4; maxiter= 100; [xO, fO] = newtons(g, xO, tolx, maxiter) 2.6. Tính đạo hàm bằng phương pháp nội suy Giả sử ta có hàm cho dưới dạng bảng X Xo Xl Xo Xn Y yo yo yo Xn Để tìm đạo hàm của hàm tại một điểm nào đó, ta nội suy hàm đó rồi tính đạo hàm cùa hàm tại điểm đã cho. Dùng hàm diffinterpO để thực hiện công việc trên, function df = diffinterp(x, y, n, xO) %Tinh dao ham cap 1 hoặc 2 bang phuong phap noi suy px = lagrange(x, y); % Tim da thuc noi suy Lagrange (x, y) [p, dp, ddp] = peval(px, xO); fprintf^Tri so cua ham la: %í\n/,p) if n ==1 df = dp; else đf = ddp; end fprintf('Dao ham cap %d la: %fW,n, df); Ví dụ 1: Để tính đạo hàm ta dùng hàm diffinterp (): clear, clc xO = pi/4; x = (2:6]*pi/16; 71 y = sin(x); X = [1.5 1.9 2.1 2.6 3.2]; y = [1.0628 1.3961 1.5432 1.8423 2.0397]; n = 2; df = diffinterp(x, y, n, xO); 2.7. Tính vi phân 2.7.1. Phương pháp chuỗi Taylor Dùng hàm TaylorO để tính vi phân: function [xout, yout] = taylor(deriv, X, y, X l , h) % Tich phan chuoi Taylor bac 4. % X, y = cac gia tri dau; i la vecto hang. % X l = gia tri cuoi cua X if size(y,l) > 1; y = y end xout = zeros(2, 1); yout = zeros(2, length(y)); xout(l) = x; yout( 1,:) = y; k = 1; while X < x l h = min(h,xl - x); d = feval(deriv,x,y); % Dao ham cua [y] hh = 1; for j = 1:4 % tao chuoi Taylor hh = hh*h/j; y = y + d(j, :)*hh; end X = X + h; k = k + 1; xout(k) = x; yout(k,:) = y; end 72 Ví dụ 1: Để giải phương trình Ta dùng hàm taylor 0 clear all, clc y = @f5; á = 0; b = 2; ya = [0 1]; h = 0.2; [x, y] = taylor(y, a, ya, b, h) plot(x, y); 2.7.2. Phuơng pháp Euler Dùng hàm euler() để thực hiện: function [X, Y] = euler(fxy, xo, xf, yo, n) %Giai phuong trinh /( x ) = f(x,y(x)) hay y’ = f(x) if n < 2 n = 2; end h = (xf - xo)/n; X = zeros(n+l, 1); M = max(size(yo));% so phuong trinh (so cot cua ma tran Y) Y = zeros(n+l, M); %dat dieu kien dau X = x o ; X (l) = x; y = yo; Y (l,:) = y/; for i = 1 :n if nargin(fxy) > 1 kl = h*feval(fxy, X, y); else kl = h*feval(fxy, x); end y = y + k l; X = X + h; X(i+1) = x; Y (i+ l,:) = y/; end function dy = f 1 (t, y) dy = zeros(3, 1); dy(l) = y(2)*y(3); dy(2) = -y (l)* y (3 ); dy(3) = -0.51 * y (l)* y (2 ); Ví dụ 1: Để giải phương trình cho bời hàm fl(x, y) ta dùng hàm euler 0 clear all, clc a = 0; b = 1; y = @ fl; ya = [0 1 1]'; m = 200; [x, y] = euler(y, a, b, ya, m) plot(x, y); 2.8. Tính tích phân 2.8.1. Công thức hình thang Dùng hàm trapezoid() để thực hiện thuật toán, function J = trapezoid(f, a, b, maxiter, tol) % Quy tac hinh thang lap. % Cu phap: J = trapezoid(f, a, b, k) fa = feval(f, a); fb = feval(f, b); J l = (fa + fb)*(b - a)/2; for k = 2:maxiter n = 2A(k -2 ); % so diem moi h = (b - a)/n; % khoang chia moi 74 X = a + h/2.0; % toa do diem moi thu nhat sum = 0.0; for i = 1 :n fx = feval(f, x); sum = sum + fx; X = X + h; end J = (J1 + h*sum)/2; if abs(Jl - J) < tol break; end J1 =J; end Ví dụ 1: Để tính tích phân ta dùng hàm trapezoid 0 clear all, clc f = inline(/(xA3+l)*sin(x)//x/); a = 0; b = 1; maxiter = 50; tol = le-6; J = trapezoid(f, a, b, maxiter, tol) 2.8.2. Công thúc Simpson Dùng hàm simpsonO để thực hiện thuật toán: function s = simpson(f, a, b, n) %n so khoang chia % neu f chua trong mot file dung ki hieu @ de goi % s = simpson(@f, a, b, n). % neu f la ham inline % s = simpson(f, a, b, n). if mođ(n, 2) ~= 0 n = n + 1 end h = (b - a)/(2*n); 75 si =0; s2 = 0; for k = 1 :n X = a + h * (2 * k -l); si = sl+ f(x); end fo rk = l:(n -l) X = a + h*2*k; s2 = s2 + f(x); end s = h*(f(a) + f(b) + 4*sl + 2*s2)/3; clc Ví dụ 1: Để tính tích phân ta dùng hàm simpson (): clear all, clc f = inline(/exp(x).*sin(x)/,V); a = 0; b = 1; n = 6; s = simpson(f, a, b, n). §3. M-FILE Trong mục trước chúng ta đã biết cách sử dụng các phép toán với các biến, để mờ rộng khả năng ứng dụng Matlab, cần sử dụng M-file nhằm hỗ trợ tương tác trong quá trình lập trình và tính toán. 3.1. Giói thiệu M-file M-file là các file ASCII (file text) chứa các (câu) lệnh Matlab, các file này có phần mờ rộng là '.m'. Có hai loại M-file: Script và Function. Các Script và Function files cũng hoạt động gần như các Procedures và Functions trong các ngôn ngữ lập trình thông dụng khác, v ề cơ bản nội dung của một script file được hiểu giống như khi nội dung đó được gõ vào tại dấu nhắc cửa sổ nhập lệnh. Hiểu đơn giản thì nó chỉ thực hiện một chuỗi các câu lệnh của Matlab. 76 ư u điểm của phương pháp sử dụng script là: - Tạo ra và xem xét, chinh sửa một chuỗi nhiều dòng lệnh (thường là 4, 5 dòng trở lên). - Có thể dễ dàng xem lại hoặc thực hiện lại công việc sau này. - Chạy các tính toán (công việc) đòi hỏi cường độ cao của CPU trên nền, xử lí kết quả, lưu lại tự động và cho phép log-off (liên quan tới UNIX). 3.2. Biên soạn và thực thi M-file + Biên soạn: Matlab cung cấp công cụ biên soạn các M-file khá tốt, đó là Matlab editor. Tuy nhiên có thể tự do sử dụng các ứng dụng soạn thảo khác cho file text của Windows như Notepad, Textpad... + Khỏi dộng Matlab Editor: Từ menu File/New/M-file, hoặc nhấn tổ hợp phím tắt ‘Ctrl - N \ hay cách nhấn vào nút ‘New Doccument’ trên thanh công cụ, cách đánh vào cửa sổ nhập lệnh ‘edit’ và tên file (nếu file chưa tồn tại trong thư mục hiện thời, Matlab sẽ hỏi để khẳng định rằng người dùng muốn tạo ra một file mới với tên như vậy). + Soạn thảo các cảu lệnh và Save: Để biết trong thư mục hiện tại (current directory) có những M-file nào, có thể sử dụng lệnh » what Để xem nội dung của một M-file, hãy nháy đúp vào file đó để mờ nó ra hoặc đánh lệnh » type tên-file + Thực thi: Để có thể thực thi một M-file, thì M-file cần phải tồn tại trong thư mục hiện thời (xem cửa sổ Current Directory), ta có thể di chuyển giữa các thư mục trong ổ cứng gẩn giống như với trình duyệt Exprorer của Windows hoặc dùng lệnh editpath (path là đường dẫn đến thư mục mà Matlab sẽ tìm kiếm file ờ đó). + Biên dịch: Không cần thiết biên dịch cả hai loại M-file của Matlab. Muốn thực hiện các lệnh chứa trong file này chỉ cẩn đánh tên file (không cần phần mờ rộng '.m') từ dấu nhắc cửa sổ lệnh. Các chỉnh sửa đã tiến hành với file và ghi lại vào ổ đĩa sẽ được thực thi khi gọi function hay script đó lần sau. Ví dụ: Gọi thực thi các lệnh có trong file baitap2.m như sau: » baitap2 Chỉ có các thông số đầu ra sẽ được thể hiện trên màn hình, chứ không phải bản thân các câu lệnh. Đế có thể xem các câu lệnh có trong file cùng lúc với các thông số dầu ra, đánh lệnh » echo on và lệnh 'echo off sẽ tắt chức nàng này. 77 3.3. Chú thích (comments) Một dạng quan trọng trong M-file là câu chú thích, được bắt đầu bằng kí tụ phần trăm (%). Bất cứ phần text nào sau kí tự '%' trên một dòng lệnh sẽ đưọc Matlab bỏ qua khồng thực hiện (trừ trường hợp kí tự % là một phẩn cùa chuỗi kí tụ giữa hai dấu nháy ‘ Mục đích chính của tính năng này là cho phép bổ sung các câu chú thích (comments) vào script file, mô tả rõ ràng mục đích, tính năng các lệnh, đoạn, vòng lặp, biến... để tiện theo dõi và dễ dàng chỉnh sửa chương trình. Ví dụ: Giả sử trong file baitap2.m có nội dung sau: % script nay giai quyet cac bai tap ve nha % lien quan toi kien thuc chuyen nganh ki thuat bien % z= rand(l); % muc nuoc bien a=omega*t*sin(2*pi); % bien do song... Khi một người dùng khác ngồi vào máy tính của bạn, muốn biết những thông tin cơ bản nhất xem file baitap2.m viết về vấn để gì, họ có thể đánh vào cửa sổ lệnh: » h e lp baitap2 và kết quả nhận dược sẽ là: ‘script nay giai quyet cac bai tap ve nha’, ‘lien quan toi kien thuc chuyen nganh ki thuat bien’. 3.4. Các hàm M-file (function M-files) Trước tiên chúng ta cần phân biệt cấc hàm M-file và các hàm số dựng sẵn, hàm trong một dòng. + Hàm dựng sẵn: sq rt(), lo g (), ex p (), sin( )... + Hàm trong một dòng (inline function): Là cách đơn giản nhất mà người dùng có thể định nghĩa một hàm. V í dụ 1: » f = inline('x*sin(x)+2’), f(5) f = Inline function: f(x) = x*sin(x)+2 ans = -2.7946 + Hàm vói M-file: Dùng cho các hàm phức tạp hơn, chẳng hạn như các vòng lặp, cáu điều kiện... bạn cẩn dùng M-file để khai báo các hàm đó. 78 Tiếp theo cẩn phân biệt các hàm M-file và các script-file: + Script M-file không phải là một hàm. Nó không có các tham số đẩu vào cũng như đầu ra và đơn giản nó chỉ thực hiện một chuỗi các câu lệnh của Matlab vói các biến được định nghĩa trong không gian làm việc. + Hàm M-file khác vói script M-file ở chỗ nó có một dòng định nghĩa hàm, qua dó liên hệ giữa các tham sô' đẩu vào và đầu ra. Hàm là cách chủ yếu dể phát huy khả năng của Matlab. So vói các script, cấc hàm có khả năng phân chia nhiệm vụ tốt hơn nhiều. Các bước chính cẩn tuân theo khi khai báo một hàm trong Matlab là: - Đặt tên cho hàm, lưu ý rằng tên đó không được xung đột với các tên đã được Matlab dành trước. Dòng đầu tiên của file này cẩn có dạng thức như sau: function[các outputs] = tên-hàm(các inputs) Một M-file hàm hoàn chỉnh có thể có dạng như sau: function [A] = đientich(a,b,c) % Tinh dien tich cua mot tam giac % khi biet chieu dai 3 canh la a, b va c. % Dau vao: % a,b,c: Chieu dai cua cac canh % Dau ra: % A: Dien tich tam giac % Cach su dung (cu phap): % Dientichcantinh = dientich(2,3,4); % Nguoi viet: Ng.Ba.Tuyen, 2007. s = (a+b+c)/2; A = sqrt(s*(s-a)*(s-b)*(s-c)); %%%%%%%%% ket thuc dientich %%%%%%%%%%% Một khíạ cạnh quan trọng khác của hàm M-file là hầu hết các hàm xây dựng trong Matlab (trừ những hàm lõi toán học) đều là các M-file mà ta có thể dọc và copy. 3.5. Đọc dữ liệu từ file và ghi ra file Nhập dữ liệu trực tiếp từ bàn phím sẽ trờ nên không thể (không thực tế) khi lượng dữ liệu quá lém. Dữ liệu đó được dùng cho phân tích nhiều lần. Trong những trường hợp này thì người sử dụng Matlab sẽ chọn cách nhập/xuất dữ liệu với file dữ liệu. Hai lệnh save và load có chức năng ghi và đọc giá trị của các biến vào/ra từ đĩa. 79 Khi làm việc với file dữ liệu, điểu cốt yếu cần lưu ý là định dạng của dữ liệu phải đúng. Định dạng dữ liệu là chìa khóa quyết định việc biên dịch dữ liệu. Có hai dạng file dữ liệu: formatted và unformatted (có định dạng và khồng định dạng). - File dữ liệu có định dạng sử dụng các định dạng chuỗi đê khai báo chính xác xem dữ liệu dược lưu ở vị trí nào. - File dữ liệu không định dạng thì khác, nó chỉ định rõ định dạng của số. V í dụ 1: Giả sử file dữ liệu dạng số được lưu với tên ‘table.dat’ trong thư mục hiện hành: 100 2256 200 4564 300 3653 400 6798 500 6432 Để sử dụng ta dùng 3 lệnh sau » fid = fopenCtable.dat','r'); » a = fscanf(fid,'%3d%4d'); » fclose(fid); Giải thích: - Mờ một file để đọc, việc này được chỉ định bằng chuỗi V , (r là viết tắt của read). Biến fid được gán cho một giá trị bằng một số nguyên tố duy nhất, đặc trưng cho file sẽ sử dụng (số này còn gọi là số chì thị của file). - Đọc vào bộ nhớ từng cặp số từ file (file có số chỉ thị là fid), một số có 3 chữ số và một số có 4 chữ số. - Đóng file (file có số chỉ thị là fid). Quá trình này tạo ra một vectơ cột chứa các phần tử 100 2256 200 4564 500 6432. 80 §4. SYMBOLIC MATH TOOLBOXE 4.1. Khái niệm chung Symbolic Math Toolboxe kết hợp tính toán bằng chữ vào môi trưìma Matlab. Các toolbox này bổ sung các tiện ích số và đồ thị với các kiểu tính loan học khác nhau. Tiện ích Nội dung Calculus Đạo hàm, tích phân, giới hạn, tống và chuỗi Taylor Linear Algebra Nghịch đảo, định thức, aiá trị riêng, phân tích và dạng chính tắc của ma trận Simplification Phương pháp rút aọn các biểu thức đại số Solution of Equations Giải bằng số, chữ các phirơns trình đại số và vi phân Variable-Precision Arithmetic Đánh giá độ chính xác cùa các biêu thức đại số Transform Biến đổi Laplace, Fourrier và z Special Mathematical Function Các hàm toán học đặc biệt cùa các ứng dụna loán học kinh điển 4.2. Khỏi động Toolbox 4.2.1. Các đối tượng ch ữ Symbolic Math Toolbox định nghĩa một kiểu dữ liệu Matlab mới aọi là đối tượng chữ hay sym. Một đối tượng chữ là một cấu trúc số liệu mù 11Ó biểu diền chuỗi các kí tự. Symbolic Math Toolbox dùng các đối tirợns chữ đê biếu diễn các biến chữ, các biểu thức chữ và các ma trận chữ. 4.2.2. Tạo các biến và các biểu thúc ch ữ Lệnh sym cho phép la xây dựng các biến và các biểu thức chữ. Ví dụ i: X = s y m ('x ’ ) a = sym('alpha') Tạo ra các biến chữ là X và a, với X là X và ;i là alpha. Ví dụ 2: Giá sir khi muốn clùns biến chữ đế hiếu diễn ti lệ: p = * + v :i , la clùnc lệnh: