- Trang Chủ
- Khoa học tự nhiên
- Tiểu luận : Kết hợp máy tính bỏ túi và Maple giải gần đúng nghiệm của bài toán Cauchy cho phương trình vi phân thường
Xem mẫu
- VIỆN TOÁN HỌC
MÔN HỌC: GIẢI TÍCH SỐ
T IỂU LUẬN
K ẾT HỢP MÁY TÍNH BỎ TÚI VÀ MAPLE GIẢI GẦN ĐÚNG NGHIỆM
CỦA BÀI TOÁN CAUCHY CHO PHƯƠNG TR ÌNH V I PHÂN THƯ ỜNG
Người thực hiện: Phạm Thị Thuỳ
Lớp: Cao học K19 - Viện Toán
HÀ NỘI – 2012
1
- Để nghiên cứu phương trình vi phân, người ta thường không giải trực tiếp
phương trình, mà sử dụng hai phương pháp: phương pháp đ ịnh tính và phương
pháp giải g ần đúng - tìm nghiệm dưới dạng xấp xỉ.
Để giải gần đúng phương trình vi phân, người ta thường dùng p hương
pháp giải tích và phương pháp số - tìm nghiệm xấp xỉ dưới d ạng các giá trị số
của nghiệm tại một số điểm trên đoạn ( a, b) và kết q uả được cho dưới dạng
bảng, như phương pháp đường gấp khúc Euler, phương pháp Runge-Kutta,...
Nhằm minh họa cho khả năng sử dụng m áy tính điện tử để giải p hương
trình vi phân, có th ể thể hiện phương pháp Euler và phương pháp Runge -Kutta
trên m áy tính điện tử khoa học Casio fx-570 ES và trên chương trình Maple qua
một số ví dụ được trình bày dưới đây.
1.1. Bà i toán Cauchy của phương trình vi phân cấp một
Một phương trình vi phân cấp một có thể viết dưới dạng giải được
y / f x, y m à ta có thể tìm được hàm y từ đạo hàm của nó. Tồn tại vô số nghiệm
tho ả mãn phương trình trên. Mỗi nghiệm phụ thuộc vào một hằng số tuỳ ý. Khi
cho trước giá trị ban đầu của y là y0 tại giá trị đầu x0 ta nhận được một nghiệm
riêng của phương trình. Bài toán Cauchy (hay bài toán có đ iều kiện đầu) tóm lại
như sau: Cho x sao cho b x a , tìm y(x) thoả mãn điều kiện
y / x f x, y
(1.1)
y x0 y0
Một cách tổng quát hơn người ta định nghĩa hệ phương trình b ậc một:
y1/ f1 x, y1 , y2 ,..., yn
/
y2 f 2 x, y1 , y2 ,..., yn
....
y / f x, y , y ,..., y
n n 1 2 n
Hệ trên có thể viết dưới dạng y / f x, y , trong đó
2
- f1 y1
f2 y2
f ... y ...
... ...
f y
n n
1.2 Giải bài toán Cauchy cho phương trình vi phân trên máy tính điện tử
và Maple
Công thức tính x ấp xỉ nghiệm theo phương pháp Euler, phương p háp
Euler cải tiến và phương pháp Runge-Kutta cho thấy, việc giải gần đúng
phương trình vi phân (1.1) có thể d ễ dàn g thực hiện tính toán trên m áy tính
khoa học Casio fx-570 ES hoặc lập trình trên Maple.
D ưới đ ây trình bày cách giải bài toán Cauchy cho một p hương trình vi
phân bằng phương pháp Euler, Euler cải tiến và p hương p háp Runge-Kutta với
các bước nội suy khác nhau trên má y tính khoa học Casio FX -570 ES và trên
Maple.
Bài 1: Sử dụng phương pháp Euler, phương pháp Euler cải tiến và phương pháp
Rungge-Kutta với độ dài bước h = 0,1 và h = 0,5 để tìm xấp xỉ nghiệm của
dy
x 2 y 2 tho ả mãn đ iều kiện ban đầu y(0) = 0 trên đo ạn 0;1 .
phương trình
dx
dy
x 2 y 2 với điều kiện ban đầu x0 =
Giải: Phải tìm nghiệm của phương trình
dx
0, y0 = 0.
Với h = 0,1 ta có:
yn 1 h. f xn , yn yn 0,1( xn yn ) yn
2 2
(1.2)
Ta có:
2 2
y1 0,1( x0 y0 ) y0 0,1(0 0) 0 0.
Với x1 = x0 + h = 0,1:
y2 0,1( x12 y12 ) y1 0.1.(0.12 0.12 ) 0 0, 001.
Tiếp tục như trên ta tính được các giá trị yn theo công thức:
2 2
yn 1 h. f xn , yn yn 0,1( xn yn ) yn .
3
- Thực hiện phép lặp (1.2) trên Casio fx -570ES:
Khai báo công thức yn1 h. f xn , yn yn 0,1( xn yn ) yn :
2 2
x2 y2
0.1 ( ÂLPH A X + ALPHA Y
) + ALPHA Y
đ ể chứa giá trị xn và dùng ô nhớ Y
Trong quy trình này, ta đã dùng ô nhớ X
để chứa giá trị của yn.
CALC để tính giá trị của yn: CALC
Dùng
Máy hỏi: X?
Khai báo: 0 và bấm phím =
Máy hỏi: Y?
Khai báo: 0 và bấm phím = (Kết quả: 0).
Kết quả trên màn hình là 0, tức là :
2 2
y1 0,1( x0 y0 ) y0 0,1(0 0) 0 0.
Đưa kết quả vào ô nhớ Y : SHIFT STO Y
Trở về công thức ban đầu (1.2): Bấm phím
Quy trình:
Tính tiếp: CALC
Máy hỏi: X?
=
Khai báo: 0.1 và bấm phím
Máy hỏi: Y? Bấm phím = (y1 = 0 vì đã sẵn có trong ô nhớ Y nên không cần
khai báo lại).
1
Kết quả trên màn hình: , tức là
1000
y2 0,1( x12 y12 ) y1 0,1(0,12 0,12 ) 0 0,13.
Y : SHIFT STO Y
Đưa kết quả vào ô nhớ
Trở về công thức ban đầu:
4
- Lặp lại quy trình với thay đổi duy nhất là khi máy hỏi X? thì khai báo các giá trị
tiếp theo: 0.2; 0.3; 0.4; …; 1.0 ta sẽ được bảng giá trị tính toán như sau:
n xn-1 yn n xn-1 yn
1 0 0 6 0,5
0,05511234067
2 0,1 0,001 7 0,6
0,09141607768
3 0,2 8 0,7
5,001 103 0,1412517676
4 0 ,3 9 0,8
0,0140026001 0,2072469738
5 10 0,9
0 ,4
0,03002220738 0,2925421046
Thực hiện phép lặp (1.2) trên Maple:
Trong Maple, đ ể tìm các giá trị yi theo công thức lặp ta có thể sử dụng mặc
định (option) remember (nhớ). Mặc định này của Maple cho phép nhớ các giá trị
cũ để tính yn, mà không cần tính lại giá trị yn-1.
Trước tiên ta khởi động chương trình Maple nhờ lệnh restart:
[> restart:
Khai báo hàm f:
[> f:=(x,y)->x^2+y^2;
f : x, y x 2 y 2
Khai báo bước nội suy h = 0,1:
[> h:=0.1;
h:=0.1
Khai báo cách tính các giá trị của x n+1 = x n + h (với x0 = 0):
[>x:=n->n*h;
x : n n h
Khai báo các giá trị ban đầu của y:
[>y(0):=0;
y(0) := 0
Khai báo thủ tục tính yn theo mặc định remember (nhớ):
5
- [>y:=proc(n) option remember;
[>y(n-1)+h*f(x(n-1),y(n-1));
[>end;
y := pro c (n ) o pti o n reme mber; y( n 1 ) hf( x( n 1 ), y( n 1 ) ) e nd
pro c
Khai báo lệnh seq (sắp xếp theo d ãy) đ ể sắp xếp các giá trị:
[>seq(y(i),i=0..10);
0, 0., 0 .00 1, 0 .00500 01 0.014002 60010 , 0 .03002220738, 0 .05511234067
, ,
0 .09141 60776 8, 0.141251767 6 , 0 .2072469738, 0 .292542 104 6
6 ,
Ta thấy kết quả này hoàn toàn trùng lặp với kết quả tính trên má y tính
khoa học Casio fx-570 ES.
Để so sánh các kết quả này với nghiệm chính xác, ta dùng lệnh dsolve (giải
phương trình vi p hân) đ ể tìm nghiệm chính xác như sau:
V ào gói công cụ Detools (công cụ Phương trìn h vi p hân):
[> with(DEtools):
Tìm nghiệm đúng của phương trình vi phân nhờ lệnh dsolve và kí h iệu n ghiệm là
Sol:
[> Sol:=dsolve({diff(Y (X),X)=X^2+Y(X)^2,Y(0)=0 },Y (X));
3 1 3 1
X BesselJ , X 2 BesselY , X 2
4 2 4 2
Sol : Y ( X )
1 1 1 1
BesselJ , X 2 BesselY , X 2
4 2 4 2
Chú ý rằng, tro ng lệnh tìm nghiệm chính xác, ta đã dùng những chữ cái in hoa để
tránh sự trùng lặp với nghiệm xấp xỉ.
Ấn định công thức nghiệm nhờ lệnh assig n:
[> a ssign(Sol);
Dùng lệnh array (lập mảng) để tạo b ảng nhằm so sánh giá trị gần đúng (tính theo
công thức Euler) và giá trị đúng của nghiệm (tính theo công thức nghiệm):
[> a rray([seq([n,y(n),evalf(subs(X=n/10,Y(X)))],n=0..10)]
6
-
0 0 0.
1
0. 0.0 0 03 3 333 4 90 60
2
0.0 0 26 6 686 9 814
0.0 0 1
3
0.0 0 90 0 347 3 190
0.0 0 50 00 1
4 0.0 1 40 0260 01 0 0.0 2 13 5 938 0 1 7
5 00 3 00 2 220 7 3 8 0.0 4 17 9 114 6 2 0
..
6 0.0 5 51 1234 06 7 0.0 7 24 4 786 1 1 8
7 0.0 9 14 1607 76 8 0.1 1 56 5985 3 6
.1 7 40 8 026 4 6
8 01 4 12 5176 7 6
.
9 0.2 0 72 4 697 3 8 0.2 5 09 0668 2 4
1 0 0.2 9 25 4 210 4 6 0.3 5 02 3 184 40
Trong b ảng nà y, cột thứ nhất là số bước lặp, các số trong cột thứ hai tương
ứng là giá trị xấp xỉ, các số trong cột thứ ba là giá trị theo công thức đúng.
Ta thấy kết quả tính toán theo công thức Euler có sai số khá lớn so với
nghiệm chính xác.
Với h = 0.05 ta có :
2 2
yn 1 h. f xn , yn yn 0, 05( xn yn ) yn
2 2
Tương tự có thể tính yn 1 h. f xn , yn yn 0,05( xn yn ) yn trên Casio
fx-570 ES bằng cách :
Khai báo công thức yn 1 h. f xn , yn yn 0,05( xn yn ) yn :
2 2
x2 y2 )
0.05 ( ÂLPH A X + ALPHA Y
+ ALPHA Y
và thao tác hoàn toàn như trên, nhưng với số bước nhiều gấp đôi (20 bước) ta
được bảng kết quả dưới đâ y.
yn yn
xn-1 xn -1
n n
1 0 0 11 0,50 0.0482462821
7
- 1
2 0 ,05 12 0,55 0.06348766728
8000
6 ,250007813
3 0 ,10 13 0,60 0.08168920148
4 0 ,15 1.750020313 10 3 14 0,65 0.1031478578
5 0 ,20 3.750173441 103 15 0,70 0.1281798318
6 0 ,25 6.875876631 103 16 0,75 0.1571263353
7 0 ,30 0.01137824052 17 0,80 0.1903607695
8 0 ,35 0.01750971373 18 0,85 0.2282976306
9 0 ,40 0.02552504324 19 0,90 0,271403621
10 0 ,45 0.03568261963 20 0,95 0.3202116173
Tính toán trên Maple:
Khai báo hàm f:
[> f:=(x,y)->x^2+y^2;
f : x, y x 2 y 2
Khai báo bước nội suy h = 0,05:
[> h:=0.05;
h:=0.05
Khai báo cách tính các giá trị của x n = x0 + n.h (với x0 = 0):
[> x:=n->n*h;
x : n n h
Khai báo các giá trị ban đầu của y:
[> y(0):=0;
y(0) := 0
Khai báo thủ tục tính giá trị yn theo công thức Euler:
[> y:=proc(n) option remember;
[> y(n-1)+h*f(x(n-1),y(n-1));
[> end;
y := pro c (n ) o pti o n remember; y( n 1 ) hf( x( n 1 ), y( n 1 ) ) e nd pro c
8
- Lập dãy giá trị của y từ 0 tới 20:
[> seq(y (i),i=1..20);
0 ., 0 .000125, 0 .000625000781,0.00175002031 ,3 ,0 .00375017344,1
0.00687587663 , 0.01137824052 0 .01750971374 0 .02552504324
1 , , ,
0 .0356826196 3 0 .04824628209 0 .0634876672 7 0 .08168920147
, , , ,
0 .103147857 8 0 .128179831 8 0 .157126335 3 0 .1903607696 0.2282976307
, , , , ,
0 .271403621 1 0 .3202116174
,
Vào gói công cụ Phương trình vi phân DEtoo ls:
[> with(DEtools):
Tìm nghiệm đúng của phương trình vi phân nhờ lệnh dsolve :
[> Sol:=dsolve({diff(Y(X),X)=X^2 +Y(X)^2,Y(0)=0},Y(X));
3 1 3 1
X BesselJ , X 2 BesselY , X 2
4 2 4 2
Sol : Y ( X )
1 1 1 1
BesselJ , X 2 BesselY , X 2
4 2 4 2
Ấn định công thức nghiệm
[> assign(Sol);
Lập mảng để so sánh giá trị gần đúng (tính theo công thức Euler) và giá trị đúng
của nghiệm (tính theo công thức nghiệm):
[> array([seq([n,y(n),evalf(subs(X=n/20,Y (X))],n=0..20 ]);
0 0
0.
1 .000041666622 1
0. 4
2
.000333334906 0
.00012 5
3 .000625000781 0 .00112502719 0
4 .00266686981 4
.00175002031 3
.0052093023 3 5
5 .00375017344 1
6 .00900347319 0
.00687587663 1
7 .0143018885 2
.0113782405 2
.0213593801 7
8 .0175097137 4
9
.0304344602 7
.0255250432 4
1 0
.0417911462 0
.035682 6196 3
1 1 .0557013376 2
.0482462820 9
.0724478611 8
1 2 .0634876672 7
1 3
.0923283103 6
.0816892014 7
1 4 .115659853 6
.103147857 8
9
- 1 5
.128179831 8 .142785233 8
1 6 .157126335 3 .174080264 6
7
1 .190360769 6 .209963219 0
1 8
.228297630 7 .250906682 4
1 9 .271403621 1 .297452631 3
0
2 .320211617 4 .350231844 0
Kết quả trùng khớp với kết quả tính toán trên Maple, có sai khác m ột đơn vị ở
chữ số thập phân thứ 10 (do làm tròn số).
Phương pháp Euler với số bứơc lặp nhiều hơn (20 bước, h = 0,05) cho kết quả
chính xác hơn;
Tính toán trên máy tính bỏ túi Casio FX-570MS bằmg phương pháp Euler
cải tiến:
Khai b áo công thức
1
f
x
xn , xn ,
y n 1 , yn h. f
h. yn f yn yn (1.3)
n 1
2
y
2
Với h = 0.1 : yn1 0,05 xn yn xn1 yn 0,1 xn yn
2 2 2 2 2
n
1
( h 0.05 và dùng lệnh CACL để tính giá trị của yn)
2
Y y2
x2 + ALPHA
0.05 ( ALPHA X + ALPHA
x2
A + Y + 0 .1 ( ALPHA X
( ALPHA
x2 y2 x2
ALPHA Y
+ ) ) ) + ALPHA Y
để chứa giá trị xn và dùng ô nhớ Y
(Trong công thức này, ta đ ã dùng ô nhớ X
để chứa giá trị của yn.
CALC
Bấm phím để tính giá trị của yn.
Máy hỏi: X?
=
Khai báo: x 0 = 0 và bấm phím
Máy hỏi: Y?
10
- Khai báo: y0 = 0 và bấm phím =
Máy hỏi: A?
=
Khai báo: 0.1 và bấm phím 0.1
1
Kết quả trên màn hình: , tức là
2000
y
2
y1 0,05 x0 y0 x12 y0 0,1 x0 y0
2 2 2 2
0
0 0, 0005.
2
0,05 02 02 0,12 0 0,1 02 02
Y: SHIFT STO Y
Đưa kết quả y1 = 0 ,0005 vào ô nhớ
Trở về công thức (1.3): Bấm phím
Tính tiếp: CALC
Máy hỏi: X?
=
Khai báo: x 1 = 0,1 và b ấm phím 0.1
Máy hỏi: Y?
Khai báo: y0 = 0 và b ấm phím = (vì y1 = 0,0005 đ ã có trong ô nhớ Y nên
không cần khai báo lại).
Máy hỏi: A?
Khai báo: 0.2 và bấm phím 0.2 =
Lặp lại quy trình với thay đổi duy nhất là khi máy hỏi X? (A?) thì khai báo các
giá trị tiếp theo: 0.1 (0.2); 0.2 (0.3); 0.3 (0.4); …; 0.9 (1.0) ta sẽ được bảng giá
trị tính toán như sau: (trùng với kết quả tính trên Maple đến chữ số cuối cùng).
xn1 yn yn
xn1
N n
1
1 0 6 0,5 0 .07344210065
200
2 0 ,1 7 0,6 0 .116816584
3.000125004 10 3
3 0 ,2 8 0,7 0 .1753963673
9.503025759 10 3
4 0 ,3 9 0,8 0 .2523742135
0.02202467595
5 0 ,4 10 0,9 0 .3518301325
0.04262140863
11
- Tính toán trên Maple :
Khởi động chương trình :
[> restart ;
Khai báo vế phải của phương trình (hàm f):
[> f:=(x,y)->x^2+y^2;
f : x, y x 2 y 2
Khai báo bước nội suy h = 0,1:
[> h:=0.1;
h:=0.1
Khai báo công thức tính xn = x0 + n.h (với x0 = 0 ):
[> x:=n->n*h;
x : n n h
Khai báo thủ tục tính giá trị yn theo công thức Euler cải tiến:
[> y:=proc(n) option remember;
[> y(n-1)+h/2*f(x(n-1),y(n-1))+f(x(n),y(n-1)+h*f(x(n-1),y(n-1))));
[> end;
y := pro c (n )
o pti o n rememb er;
y( n 1 )
1/2h ( f( x( n 1 ), y( n 1 ) ) f( x( n ), y( n 1 ) h f( x( n 1 ), y( n 1 ) ) ) )
e nd pro c
Khai báo giá trị ban đầu của y:
[> y(0):=0;
y(0) := 0
Lập dãy giá trị của y từ 0 tới 10:
[> seq(y (i),i=0..10);
0, .000500000000,0.00300012500, .00950302575, .02202467594
4 9 ,
.0426214086 3.0734421006 5 .1168165840 .175396367 3 .252374213 4
, , , , ,
.351830132 5
Vào gói công cụ Phương trình vi phân DEtoo ls:
12
- [> with(DEtools):
Tìm nghiệm đúng của p hương trình vi phân nhờ lệnh dsolve:
[> Sol:=dsolve({diff(Z(X),X)= X^2+(Z(X))^2,Z(0)=0},Z(X));
3 1 3 1
X BesselJ , X 2 BesselY , X 2
4 2 4 2
Sol : Z ( X )
1 1 1 1
BesselJ , X 2 BesselY , X 2
4 2 4 2
Ấn định công thức nghiệm
[> assign(Sol);
Lập mảng để so sánh giá trị gần đúng (tính theo công thức Euler) và giá trị đúng
của phương trình (tính theo công thức nghiệm ):
[> array([seq([n,y(n),evalf(subs(X=n/10,Z(X)))],n=0..10]);
0 0 0.
1
.0 0 0500 0 0000 0 0 .0 0 0333 3 3490 6 0
2
.0 0 3000 1 2500 4 .0 0 2666 8 6981 4
3
.0 0 9503 0 2575 9 .0 0 9003 4 7319 0
4 .0 2 1359 3 801 7
.0 2 2024 6 759 4
.0 4 1791 1 462 0
5 .0 4 2621 4 086 3
6 .0 7 2447 8 611 8
.0 7 3442 1 006 5
7 .1 1 5659 8 53 6
.1 1 6816 5 84 0
.1 7 4080 2 64 6
8 .1 7 5396 3 67 3
9 .2 5 0906 6 82 4
.2 5 2374 2 13 4
0
1 .3 5 1830 1 32 5 .3 5 0231 8 44 0
Kết q uả tính toán trê n Casio fx-570 ES hoàn toàn trùng khớp với kết quả
tính toán trên Maple. H ơn nữa, chỉ cần với h=0.1, p hương pháp Euler cải tiến đã
cho kết quả tốt hơn phương p háp Euler với h=0.05.
Tương tự, ta cũng đi tính xấp xỉ nghiệm nhờ phương pháp Euler cải tiến trên
Map le khi h=0,05 như sau.
Khởi động chương trình:
[> restart;
Khai báo vế phải của phương trình (hàm f ):
13
- [> f:=(x,y)->x^2+y^2;
f : x, y x 2 y 2
Khai báo bước nội suy h = 0,05:
[> h:=0.05;
h:=0.05
Khai báo công thức tính xn = x0 + n.h (với x0 = 0 ):
[> x:=n->n*h;
x : n n h
Khai báo thủ tục tính giá trị yn theo công thức Euler cải tiến:
[> y:=proc(n) option remember;
[> y(n-1)+h/2*f(x(n-1),y(n-1))+f(x(n),y(n-1)+h*f(x(n-1),y(n-1))));
[> end;
y := pro c (n )
o pti o n rememb er;
y( n 1 )
1/2h ( f( x( n 1 ), y( n 1 ) ) f( x( n ), y( n 1 ) h f( x( n 1 ), y( n 1 ) ) ) )
e nd pro c
Khai báo giá trị ban đầu của y:
[> y(0):=0;
y(0) := 0
Lập dãy giá trị của y từ 0 tới 20:
[> seq(y (i),i=0..20);
0, .0000625000 000,0.0003750009 76,8.00118752363,4.00275019259 ,
2
.00531344588 , .00912843247 , .01444766188 .0215259718 5 .0306218 8483
0 8 , , ,
.04199943062.05593052466 .07269800874 .0925994870 6 .1159521276
, , , , ,
.143098652 2.174414813 0 .210318759 0 .2512828469 .2978486637
, , , , ,
.350646340 8
Vào gói công cụ Phương trình vi phân DEtoo ls:
[> with(DEtools):
Tìm nghiệm đúng của p hương trình vi phân nhờ lệnh dsolve:
[> Sol:=dsolve({diff(Z(X),X)= X^2+(Z(X))^2,Z(0)=0},Z(X));
14
- 3 1 3 1
X BesselJ , X 2 BesselY , X 2
4 2 4 2
Sol : Z ( X )
1 1 1 1
BesselJ , X 2 BesselY , X 2
4 2 4 2
Ấn định công thức nghiệm
[> assign(Sol);
Lập mảng để so sánh giá trị gần đúng (tính theo công thức Euler) và giá trị đúng
của phương trình (tính theo công thức nghiệm ):
[> array([seq([n,y(n),evalf(subs(X=n/20,Z(X)))],n=0..20]);
0 0 0.
.0000625000000 0 .000041666622 1
1 4
.000375000976 8 .000333334906 0
2
.00112502719 0
3 .00118752363 4
.00266686981 4
4 .00275019259 2
.0052093023 3 5
5 .00531344588 0
6
.00900347319 0
.00912843247 8
7 .0143018885 2
.0144476618 8
8
.0215259718 50
.0306218848 3 7 .0304344602 7
0213593801
9
.
1 0
.0417911462 0
.0419994306 2
1 1 .0557013376 2
.0559305246 6
.0724478611 8
1 2 .0726980087 4
1 3 .0923283103 6
.0925994870 6
1 4 .115659853 6
.115952127 6
1 5 .142785233 8
.143098652 2
1 6
.174080264 6
.174414813 0
1 7 .209963219 0
.210318759 0
1 8 .2 50906682 4
.251282846 9
1 9 .297452631 3
.297848663 7
2 0
.350646340 8 .350231844 0
Kết quả tính toán trên Casio fx-5 70 ES hoàn toàn trùng khớp với kết quả tính
toán trên Maple. V ới cùng số bước lặp (n=20, h=0.05), phương pháp Euler cải
tiến cho kết quả tốt hơn phương pháp Euler rất nhiều.
Phương phá p Rung e-Kutta cấp bốn
Ta có: f(x,y) = x2 +y2, x 0 = 0, y0 = 0, áp dụng công thức ta được :
15
- 2 2
k1 f xn , yn xn yn
2 2
h hk 0.1 0.1k1
k 2 f xn , yn 1 xn yn
2 2 2 2
2 2
h hk 0.1 0.1k2
k3 f x n , y n 2 xn yn
2 2 2 2
2
2
k4 f xn1 , yn hk3 xn 1 yn 0.1k3
h 0.1
k1 2k2 2k3 k4 yn k1 2k 2 2k3 k 4
và yn 1 yn
6 6
Khởi động chương trình:
[> restart ;
Định nghĩa yrk (tính y theo Runge-Kutta):
[> yrk:='yrk';
yr k : = y r k
Khai báo vế phải của phương trình (hàm f ):
[> f:=(x,y)->x^2+y^2;
f : x, y x 2 y 2
Khai báo bước nội suy h = 0,1:
[> h:=0.1;
h:=0.1
Khai báo công thức tính xn = x0 + n.h (với x0 = 0 ):
[> x:=n->n*h;
x : n n h
Khai báo thủ tục tính giá trị yn theo công thức Runge-Kutta cấp bốn:
[> yrk:=proc(n)
[> local k1,k2,k3,k4;
[> option rememb er;
[> k1:=f(x(n-1),yrk(n-1));
[> k2:=f(x(n-1)+h/2,yrk(n-1)+h*k1/2);
[> k3:=f(x(n-1)+h/2,yrk(n-1)+h*k2/2);
[>k4:=f(x(n),yrk(n-1)+h* k3);
[> yrk(n-1)+h/6*(k1+2*k2+2*k3+k4)
[> end;
yrk := pro c (n )
16
- l o c al k1, k2, k3, k4;
o pti o n remember;
k1 := f( x( n 1 ), yrk( n 1 ) ) ;
k2 := f( x( n 1 ) 1/2h , yrk( n 1 ) 1/2h k1 ) ;
k3 := f( x( n 1 ) 1/2h , yrk( n 1 ) 1/2h k2 ) ;
k4 := f( x( n ), yrk( n 1 ) hk3 ) ;
yrk( n 1 ) 1/6h ( k 1 2k2 2k3 k4 )
e nd pro c
Khai báo giá trị ban đầu của y:
[> y(0):=0;
y(0) := 0
Lập dãy giá trị của y từ 0 tới 10:
[> seq(yrk(i),i=0..10);
0, .000333334895,8.00266687536, .0090034981 3, .02135944733
9 1 ,
.0417912884 8 .07244812485 .1156603048 .1740810040 .2509078684
, , , , ,
.350233741 7
Vào gói công cụ Phương trình vi phân DEtoo ls:
[> with(DEtools):
Tìm nghiệm đúng của p hương trình vi phân nhờ lệnh dsolve:
[> Sol:=dsolve({diff(Z(X),X)= X^2+(Z(X))^2,Z(0)=0},Z(X));
3 1 3 1
X BesselJ , X 2 BesselY , X 2
4 2 4 2
Sol : Z ( X )
1 1 1 1
BesselJ , X 2 BesselY , X 2
4 2 4 2
Ấn định công thức nghiệm
[> assign(Sol);
Lập mảng để so sánh giá trị gần đúng (tính theo công thức Euler) và giá trị đúng
của phương trình (tính theo công thức nghiệm ):
[> array([seq([n,y(n),evalf(subs(X=n/10,Z(X)))],n=0..10]);
0 0 0.
1 .0 0 0333 3 3489 5 8 .0 0 0333 3 3490 6
0
2
.0 0 2666 8 7536 9 .0 0 2666 8 6981 4
3
.0 0 9003 4 9813 1 .0 0 9003 4 7319 0
4
.0 2 1359 4 473 3 .0 2 1359 3 801 7
5 .0 4 1791 2 884 8 .0 4 1791 1 462 0
17
-
6 .0 7 2447 8 611 8
.0 7 2448 1 248 5
7
.1 1 566 0 3 04 8 .1 1 5659 8 53 6
.1 7 4080 2 64 6
8 .1 7 4081 0 04 0
9 .2 5 0906 6 82 4
.2 5 0907 8 68 4
0
1 .3 5 0233 7 41 7 .3 5 0231 8 44 0
So sánh các kết quả của phương pháp Runge-Kutta cấp 4 trong bảng trên
với kết quả đã thực hiện theo phương pháp Euler và phương pháp Euler cải
tiến, ta thấy rằng phương pháp nà y cho kết quả chính x ác hơn tại mỗi điểm so
với phương pháp Euler và phương pháp Euler cải tiến. Với số bước ít (n=10,
h=0.1) ta đã thu được kết quả tốt hơn phương pháp Euler cải tiến với số bước
gấp đôi (n=20, h=0.05).
H oàn toàn tương tự (với tha y đổi duy nhất trong chương trình là khai báo
lại b ước nội suy h=0.05), ta có thể tính theo phương pháp Runge -Kutta với số
bước n=20 (h=0.05) như sau.
Khởi động chương trình:
[> restart;
Định nghĩa yrk ( tính y theo Runge-Kutta):
[> yrk:='yrk';
y r k : = y rk
Khai báo vế phải của phương trình (hàm f ):
[> f:=(x,y)->x^2+y^2;
f : x, y x 2 y 2
Khai báo bước nội suy h = 0,05:
[> h:=0.05;
h:=0.05
Khai báo công thức tính xn = x0 + n.h (với x0 = 0 ):
[> x:=n->n*h;
x : n n h
Khai báo thủ tục tính giá trị yn theo công thức Runge-Kutta cấp bốn:
18
- [> yrk:=proc(n)
[> local k1,k2,k3,k4;
[> option rememb er;
[> k1:=f(x(n-1),yrk(n-1));
[> k2:=f(x(n-1)+h/2,yrk(n-1)+h*k1/2);
[> k3:=f(x(n-1)+h/2,yrk(n-1)+h*k2/2);
[>k4:=f(x(n),yrk(n-1)+h* k3);
[> yrk(n-1)+h/6*(k1+2*k2+2*k3+k4)
[> end;
yrk := pro c (n )
l o c al k1, k2, k3, k4;
o pti o n remember;
k1 := f( x( n 1 ), yrk( n 1 ) ) ;
k2 := f( x( n 1 ) 1/2h , yrk( n 1 ) 1/2h k1 ) ;
k3 := f( x( n 1 ) 1/2h , yrk( n 1 ) 1/2h k2 ) ;
k4 := f( x( n ), yrk( n 1 ) hk3 ) ;
yrk( n 1 ) 1/6h ( k 1 2k2 2k3 k4 )
e nd pro c
Khai báo giá trị ban đầu của y:
[> y(0):=0;
y(0) := 0
Lập dãy giá trị của y từ 0 tới 20:
[> seq(yrk(i),i=0..20);
0, .00004166667 88,7.0003333349 6 3,7.00112502731,6.00266687038 ,
2
.00520930346 , .00900347509,2.01430189176 .02135938501 .03043446755
2 , , ,
.04179115619 .05570135121 .07244787939 .09232833422 .1156598841
, , , , ,
.142785273 2 .174080314 6 .2099632826 .250906762 3 .2974527325
, , , , ,
.350231972 4
Vào gói công cụ Phương trình vi phân DEtoo ls:
[> with(DEtools):
Tìm nghiệm đúng của p hương trình vi phân nhờ lệnh dsolve:
[> Sol:=dsolve({diff(Z(X),X)= X^2+(Z(X))^2,Z(0)=0},Z(X));
3 1 3 1
X BesselJ , X 2 BesselY , X 2
4 2 4 2
Sol : Z ( X )
1 1 1 1
BesselJ , X 2 BesselY , X 2
4 2 4 2
19
- Ấn định công thức nghiệm
[> assign(Sol);
Lập mảng để so sánh giá trị gần đúng (tính theo công thức Euler) và giá trị đúng
của phương trình (tính theo công thức nghiệm ):
[> array([seq([n,y(n),evalf(subs(X=n/20,Z(X)))],n=0..20]);
0 0 0.
1 .0000625000000 0
.000041666622 1 4
2 .000375000976 8 .00 03333349060
3
.00112502719 0
.00118752363 4
4
.00266686981 4
.00275019259 2
5 .0052093023 35
.00531344588 0
6
.00900347319 0
.00912843247 8
7 .0143018885 2
.0144476618 8
.0213593801 7
8 .0215259718 5
9 .0304344602 7
.0306218848 3
1 0 .0417911462 0
.0419994306 2
1 1
.0559305246 6 .0557013376 2
1 2 .0726980087 4 .0724478611 8
1 3 .0923283103 6
.0925994870 6
1 4 .115659853 6
.115952127 6
1 5
.143098652 2 .142785233 8
.174080264 6
1 6 .174414813 0
1 7 .209963219 0
.210318759 0
1 8
.251282846 9 .250906682 4
1 9 .297848663 7 .297452631 3
2 0 .350231844 0
.350646340 8
Các kết quả của phương pháp Runge-Kutta cấp 4 là tốt hơn rất nhiều so
với kết quả đã thực hiện theo phương pháp Euler và p hương p háp Euler cải tiến
với cùng số b ước (n=20, h=0.05) và tốt hơn phương pháp Runge -Kutta với
số bước ít hơn (n=10, h=0.1).
Các thao tác này có thể được coi là các chương trình mẫu để giải các bài
to án khác (chỉ cần khai báo lại phương trình cần giải).
Bài 2: Sử dụng phương pháp Euler và phương pháp Euler cải tiến với độ dài
20
nguon tai.lieu . vn