Giải hệ phương trình phi tuyến bằng MATLAB 04/203 Nguyễn Đ.M. Anh. Giới thiệu về lệnh fsolve ở MATLAB Lệnh "fsolve" ở MATLAB là một cách nhanh chóng và hiệu quả để tìm một nghiệm của hệ phương trình phi tuyến (n phương trình và n ẩn số). Vì phương trình tuyến tính là trường hợp đặc biệt của phương trình phi tuyết nên lệnh fsolve có thể được sử dụng để tìm nghiệm của hệ phương trình tuyến tính. Bản chất của fsolve là sử dụng phương pháp số để tìm nghiệm. 2. Ví dụ về quy trình sử dụng lệnh fsolve Xem xét ví dụ sau: Đây là hệ phương trình phi tuyến tính với 2 ẩn x và y. 2. Quy trình: Phương trình này có thể được giải bằng lệnh fsolve như sau: Bước : Tạo ra một m.file (MATLAB file) để mô tả về hàm: function F = myfun(x) F = [4*x() - x(2) - exp(-2*x()); -2*x() + 2*x(2) - exp(-x(2))]; File này được lưu dưới tên: myfun.m (Chú ý add to path file này, nếu không MATLAB sẽ không nhận dạng được hàm myfun) Bước 2: Sử dùng lệnh fsolve để giải hệ phương trình. Để sử dụng lệnh fsolve, cần thực hiện 2 thao tác: + Trước hết, cần phải đoán giá trị ban đầu cho nghiệm x và x 2. Ví dụ: x 0 = - và x 20 = - (Xem thảo luận ở phần 3 về giá trị ban đầu) + Sau đó, fsolve sẽ sử dụng thuật toán để tìm kiếm nghiệm, xuất phát từ giá trị được đề xuất ban đầu. x0 = [-; -]; % Giá trị ban đầu [x,fval,exitflag] = fsolve(@myfun,x0) 2.2 Kết quả như sau: 0.2836 0.567.0e-08 * -0.99 Xem chi tiết thuật toán được sủ dụng tại: http://www.mathworks.co.uk/help/optim/ug/fsolve.html -0.0826 2.3 Chú thích về kết quả: 2 x: nghiệm của hệ phương trình (trong ví dụ trên, x là một véc tơ của 2 biến x và x 2, do vậy kết quả x sẽ được viết dưới dạng một véc tơ). Trong đó x = 0.2836 và x 2 = 0.567 fval: giá trị của hàm F (trong ví dụ trên F cũng là một véc tơ). Khi x=[0.2836; 0.567] thì F xấp xỉ 0. exitflag: giải thích tại sao thuật toán dừng (Kết quả mong muốn là vì khi này hàm hội tụ tại một nghiệm x) MATLAB cũng cung cấp một số lựa chọn tùy theo mục đích của người nghiên cứu. Khi ấy, cần viết thêm optimset vào lệnh trên. Ví dụ: x0 = [-; -]; [x,fval,exitflag] = fsolve(@myfun,x0,optimset('display','iter')) Chi tiết về các lệnh có thể xem tại website ở ghi chú. 3. Lưu ý về giá trị ban đầu: Giá trị ban đầu càng gần với nghiệm thì thời gian để thuật toán tìm càng nhanh. Trong ví dụ trên, giá trị ban đầu của x được chọn là [-; -]. Tuy nhiên, người lập 2 Xem thêm giải thích về kết quả tại địa chỉ website ở ghi chú. trình có thể chọn các giá trị tùy ý khác như x=[0; 0] hay x=[; ], kể cả x=[-20;-20] cũng sẽ tìm ra được nghiệm x = 0.2836 và x 2 = 0.567 và. Nhưng nếu chọn x=[-50; -50], có thể sẽ gặp vấn đề sau: fsolve stopped because it exceeded the function evaluation limit, options.maxfunevals = 200 (the default value). -7.0000-7.090.0e+4 * -5.8348-0.0000 0 Trước hết, exitflag=0, nghĩa là hàm chưa hội tụ tại nghiệm tìm được (Dễ thấy fval khác với 0). Vấn đề này do giới hạn định giá hàm được lập định sẵn là 200, do vậy quy trình tìm kiếm nghiệm dừng lại sau 200 lần. Trong tình huống này, có hai cách xử lý: - Thay đổi giá trị ban đầu x0. - Thay đổi giới hạn đánh giá hàm từ 200 đến, ví dụ như, 0000. Để tiến hành điều này, sử dụng optimset như đã để cập ở trên. Ví dụ như: x0 = [-50; -50]; [x,fval,exitflag] = fsolve(@myfun,x0,optimset('display','iter', 'MaxFunEvals', 0000 )) Kết quả thu được tương tự như kết quả đã đạt được trên: 0.2836 0.567.0e- * -0.373-0.46 4. Sử dụng lệnh fsolve để giải hệ phương trình tuyến tính Ví dụ: Viết lại hệ phương trình như sau: Bước : Viết hàm fun.m function F = fun(x) F = [2*x() + 3*x(2) + x(3)-7; x() + 4*x(2) - 2*x(3)-; x() - x(2) + x(3)-2]; Bước 2: Tìm nghiệm x0 = [0; 0; 0]; % Giá trị ban đầu [x,fval,exitflag] = fsolve(@fun,x0) % Call solver Kết quả: 2 0 0 0 5. Vấn đề nghiệm duy nhất và đa nghiệm Một điều quan trọng nên lưu ý khi sử dụng lệnh fsolve là: Lệnh fsolve chỉ cung cấp một nghiệm của hệ phương trình, dựa trên giá trị ban đầu được cung cấp. Ví dụ: Xem xét hệ phương trình trên, ta thấy rằng phương trình thứ hai chính là kết quả của việc nhân phương trình thứ nhất với 2. Do vậy, hai phương trình này đồng nhất. function F = myfun(x) F = [4*x() - x(2) - exp(-2*x()); 8*x() - 2*x(2) - 2*exp(-2*x())]; x0 = [-; -]; % Giá trị ban đầu [x,fval,exitflag] = fsolve(@myfun,x0) Nếu x0=[-;-], kết quả thu được là: -0.035 -.087 .0e- * -0.4593-0.985 Nghiệm này phù hợp với hệ phương trình Tuy nhiên, nếu x0=[0,0] x0 = [0; 0]; % Giá trị ban đầu [x,fval,exitflag] = fsolve(@myfun,x0) Kết quả thu được là: 2.948.7645.0e-08 * -0.2376-0.4753 Nghiệm này cũng phù hợp với hệ phương trình Do vậy, khi sử dụng fsolve, nên lưu ý khả năng hệ có đa nghiệm. Để kiểm tra xem nghiệm có phải nghiệm duy nhất không, người lập trình nên sử dụng các giá trị ban đầu khác nhau để xem nghiệm tìm được có như nhau không. |
