Giáo trình Matlab cơ bản - Chương 7: Các phương trình vi phân thường

; 6 
    n = 100;  
end 
dx0(1) = (xf ‐ x0)/(tf ‐ t0); % cho gia tri dau cua xʹ(t0) 
y0 = [x0 dx0(1)]ʹ; 
[t, x] = rungekutta(f, t0, tf , y0, n); % khoi gan bg RK4 
 387
e(1) = x(end,1) ‐ xf;  
dx0(2) = dx0(1) ‐ 0.1*sign(e(1)); 
for k = 2: kmax‐1 
    y1 = [x0 dx0(k)]ʹ; 
    [t, x] = rungekutta(f, t0, tf, y1, n); 
    %sai so giua gia tri cuoi va dich 
    e(k) = x(end, 1) ‐ xf; % x(tf)‐ xf 
    ddx = dx0(k) ‐ dx0(k ‐ 1);  
    if abs(e(k))< tol | abs(ddx)< tol 
        break;  
    end 
    deddx = (e(k) ‐ e(k ‐ 1))/ddx;  
    dx0(k + 1) = dx0(k) ‐ e(k)/deddx;  
end 
Để giải phương trình: 
  ′′ ′= +2y 2y 4xyy  với điều kiện biên: y(0) = 0.25, y(1) = 1/3  
Đặt:  ′ = 1y y ,  ′′ = 2y y  ta đưa phương trình về hệ phương trình vi phân cấp 1: 
=⎧⎨ = +⎩
1 2
2 1 2 1
y y
y 2y 4xy y
và biểu diễn nó bằng hàm f7(): 
function dx = f7(t, x) %Eq.(6.6.5) 
dx(1) = x(2);  
dx(2) = (2*x(1) + 4*t*x(2))*x(1); 
Để giải phương trình ta dùng chương trình ctbvp2shoot.m: 
clear all, clc 
t0 = 0;  
tf = 1;  
x0 = 1/4;  
xf = 1/3; % thoi gian dau/cuoi va cac vi tri 
n = 100;  
tol = 1e‐8;  
 388
kmax = 10; 
y = @f7; 
[t, x] = bvp2shoot(y, t0, tf, x0, xf, n, tol, kmax); 
xo = 1./(4 ‐ t.*t);  
err = norm(x(:,1) ‐ xo)/(n + 1) 
plot(t,x(:, 1))  
3. Phương pháp sai phân hữu hạn: Ta xét phương trình: 
′′ ′= = α = βy f(x,y,y ) y(a) ,y(b)  
Ý tưởng của phương pháp này là chia đoạn [a, b] thành n đoạn nhỏ có bước h 
và xấp xỉ các đạo hàm bằng các sai phân: 
  + −′ = i 1 iy yy
2h
  − +− +′′ = i 1 i i 12y 2y yy h  
Như vậy: 
  − + +− + −⎛ ⎞′′ = = ⎜ ⎟⎝ ⎠
i 1 i i 1 i 1 i
i i2
y 2y y y yy f x ,y ,
h 2h
  i = 1, 2, 3,... 
Phương  trình bnày sec đưa đến hệ phương  trình đối với xi, yi. Ta xây dựng 
hàm bvp2fdf(): 
function [t, x] = bvp2fdf(a1, a0, u, t0, tf, x0, xf, n) 
% Giai pt : xʺ + a1*x’ + a0*x = u with x(t0) = x0, x(tf) = xf 
% bang pp sai phan huu han 
h = (tf ‐ t0)/n;  
h2 = 2*h*h; 
t = t0 + [0:n]ʹ*h; 
if ~isnumeric(a1) 
    a1 = a1(t(2: n));  
else 
    length(a1) == 1 
    a1 = a1*ones(n ‐ 1, 1); 
end 
if ~isnumeric(a0) 
    a0 = a0(t(2:n));  
else length(a0) == 1 
 389
    a0 = a0*ones(n ‐ 1, 1); 
end 
if ~isnumeric(u) 
    u = u(t(2:n));  
elseif length(u) == 1 
    u = u*ones(n‐1,1); 
else  
    u = u(:); 
end 
A = zeros(n ‐ 1, n ‐ 1);  
b = h2*u; 
ha = h*a1(1);  
A(1, 1:2) = [‐4 + h2*a0(1) 2 + ha]; 
b(1) = b(1) + (ha ‐ 2)*x0; 
for m = 2:n ‐ 2  
    ha = h*a1(m);  
    A(m, m ‐ 1:m + 1) = [2‐ha   ‐4+h2*a0(m)   2+ha]; 
end 
ha = h*a1(n ‐ 1);  
A(n ‐ 1, n ‐ 2:n ‐ 1) = [2‐ha   ‐4+h2*a0(n ‐ 1)]; 
b(n ‐ 1) = b(n‐1) ‐ (ha+2)*xf; 
x = [x0 trid(A,b)ʹ xf]ʹ; 
  function x = trid(A, b) 
% giai he pt trdiagonal 
n = size(A, 2); 
for m = 2:n  
    tmp = A(m, m ‐ 1)/A(m ‐ 1, m ‐ 1); 
    A(m, m) = A(m, m) ‐A(m ‐ 1, m)*tmp; 
    A(m, m ‐ 1) = 0; 
    b(m, :) = b(m, :) ‐b(m ‐ 1, :)*tmp; 
end 
x(n, :) = b(n, :)/A(n, n); 
for m = n ‐ 1: ‐1: 1  
    x(m, :) = (b(m, :) ‐A(m, m + 1)*x(m + 1))/A(m, m); 
end 
 390
Để giải phương trình: 
  ′′ ′+ − =22 2y y 0x x  
với y(1) = 5 và y(2) = 3 ta dùng chương trình ctbvp2fdf.m:  
clear, clc 
x0 = 1;%toa do bien dau  
y0 = 5; %gia tri bien dau 
xf = 2; %toa do bien cuoi 
yf = 3; %gia tri bien cuoi 
n = 100;%so buoc tinh 
a1 = inline(ʹ2./xʹ, ʹxʹ);  
a0 = inline(ʹ‐2./x./xʹ,ʹxʹ);  
u = 0;%ve phai cua phupng trinh  
[x, y] = bvp2fdf(a1, a0, u, x0, xf, y0, yf, n); 
plot(x, y)  
§11. PHƯƠNG PHÁP LẶP PICARD 
  Việc giải bài toán Cauchy: 
  ′ =y f(t,y) ,  =o oy(t ) y                 (1) 
hoàn toàn tương đương với việc tìm nghiệm y(t) của phương trình tích phân: 
  [ ]= + ∫1
o
t
o
t
y(t) y f z,y(z) dz               (2) 
Nội dung  của phương  trình Picard  là  thay  cho  việc  tìm  nghiệm  đúng  của 
phương trình (2) ta tìm nghiệm gần đúng theo công thức: 
  [ ]+ = + ∫1
o
t
n 1 o
t
y y f z,y(z) dz               (3) 
Trong đó để đơn giản ta chọn nghiệm gần đúng đầu tiên là: 
  yo(x) = yo 
Ta xây dựng hàm picard() để thực hiện thuật toán trên: 
function g = picard(f, y0, maxiter) 
syms x y 
g = subs(f, y0); 
 391
g = int(g, x); 
for i = 1:maxiter 
    g = subs(f, g); 
    g = int(g, x); 
end 
g = sort(g);  
Để  giải  phương  trình  ta  dùng  chương  trình 
ctpicard.m: 
clear all, clc 
syms x y 
%f = 1 + y^2;%Phuong trinh yʹ = 1+y^2 
%y0 = 0;%dieu kien dau 
f = y ‐ x; 
y0 = 2; 
g = picard(f, y0, 4) 
§12. PHƯƠNG PHÁP GILL 
  Phương pháp Gill cũng tương tự như phương pháp Runge ‐ Kutta, giá 
trị nghiệm tại yn+1 được tính bằng: 
  +
⎡ ⎤= + + − + + +⎢ ⎥⎦⎣n 1 n 1 2 4 4
1y y k (2 2)k (2 2)k k
6
Với: 
  =1 n nk hf(x ,y )  
       ⎡ ⎤= + +⎢ ⎥⎣ ⎦2 n n 1
1 1k hf x h,y k
2 2
  ( )⎡ ⎤⎛ ⎞= + + − + + −⎢ ⎥⎜ ⎟⎝ ⎠⎣ ⎦3 n n 1 21 1 2k hf x h,y 1 2 k 1 k2 2 2  
⎡ ⎤⎛ ⎞= + − + +⎢ ⎥⎜ ⎟⎝ ⎠⎣ ⎦4 n n 2 3
2 2k hf x h,y k 1 k
2 2
Ta xây dựng hàm gill() để thực hiên thuật toán trên: 
function [x, y] = gill(f, a, b, y0, n) 
%Phuong phap Gill de giai phuong trinh yʹ(x) = f(x,y(x)) 
if nargin < 4 | n <= 0 
 392
    n = 100;  
end 
if nargin < 3 
    y0 = 0;  
end 
y(1,:) = y0(:)ʹ;  
h = (b ‐ a)/n;  
x = a + [0:n]ʹ*h; 
if nargin(f) >1 
    for k = 1:n 
        k1 = h*feval(f, x(k), y(k, :));  
        k1 = k1(:)ʹ;  
        k2 = h*feval(f, x(k)+h/2, y(k, :)+k1/2);  
        k2 = k2(:)ʹ;        
        k3 = h*feval(f, x(k)+h/2, y(k, :)+(sqrt(2)‐1)*k1/2+(1‐sqrt(2)/2)*k2);  
        k3 = k3(:)ʹ; 
        k4 = h*feval(f, x(k)+h, y(k, :) ‐sqrt(2)/2*k2+(1+sqrt(2)/2)*k3);  
        k4 = k4(:)ʹ;  
        y(k+1, :) = y(k, :) + (k1 + (2 ‐ sqrt(2))*k2 + (2+sqrt(2))*k3 + k4)/6;  
    end 
else 
    for k = 1:n 
        f1 = h*feval(f, x(k));  
        f1 = f1(:)ʹ;  
        f2 = h*feval(f, x(k) + h/2);  
        f2 = f2(:)ʹ;        
        f3 = h*feval(f, x(k) + h/2);  
        f3 = f3(:)ʹ; 
        f4 = h*feval(f, x(k) + h);  
        f4 = f4(:)ʹ;  
        y(k+1, :) = y(k, :) + (f1 + (2 ‐ sqrt(2))*f2 + (2+sqrt(2))*f3 + f4)/6;  
    end 
end 
Để giải phương trình ta dùng chương trình ctgill.m: 
 393
clear all, clc 
a = 0; 
b = 1; 
y = inline(ʹx + yʹ); 
ya = 0.5; 
n = 10;%so lan tinh chi n = 10 
[t, u] = gill(y, a, b, ya, n); 
plot(t, u); 
[l, v] = rungekutta(y, a, b, ya, n); 
hold on 
plot(l, v, ʹ.rʹ) 
§13. PHƯƠNG PHÁP RUNGE – KUTTA – FEHLBERG 
  Một   phương pháp để bảo đảm độ chính xác của nghiệm của phương 
trình vi phân là giả bài toán hai lần với bươc tính là h và 0.5h rồi so sánh kết 
quả tại các nút  ứng với h. Tuy nhiên điều này đòi hỏi tăng số lần tính và phải 
tính  lại khi giá trị tại các nút khác nhau nhiều. Dùng phương pháp Runge – 
Kutta – Fehlberg có thể tránh được khó khăn này. Nó có thủ tục để xác định 
liệu kích thước bước tính h đã thích hợp chưa. Tại mỗi bươc tính, hai xấp xỉ 
khác nhau của nghiệm được  tính và so sánh với nhau. Nếu chúng sai khác 
trong phạm vi sai số cho trước thì nghiệm được chấp nhận. Nếu sai khác còn 
lơn hơn sai số cho phép thì bước h được giảm và ta lại tìm nghiệm với h mới. 
Mỗi bước tính theo phương pháp Runge – Kutta – Fehlberg đòi hỏi 6 giá trị: 
  1 i ik hf(t ,y )=   
  2 i i 1
1 1k hf t h,y k
4 4
⎛ ⎞= + +⎜ ⎟⎝ ⎠  
  3 i i 1 2
3 3 9k hf t h,y k k
8 32 32
⎛ ⎞= + + +⎜ ⎟⎝ ⎠  
  4 i i 1 2 3
12 1932 7200 7296k hf t h,y k k k
13 2197 2197 2197
⎛ ⎞= + + − +⎜ ⎟⎝ ⎠  
  5 i i 1 2 3 4
439 3680 845k hf t h,y k 8k k k
216 513 4104
⎛ ⎞= + + − + −⎜ ⎟⎝ ⎠  
  5 i i 1 2 3 4 5
1 8 3544 1859 11k hf t h,y k 2k k k k
2 27 2565 4104 40
⎛ ⎞= + − + − + −⎜ ⎟⎝ ⎠  
Xấp xỉ nghiệm theo phương pháp Runge – Kutta bậc 4: 
 394
i 1 i 1 3 4 5
25 1408 2197 1y y k k k k
216 2565 4104 5+
= + + + −  
Và nghiệm tốt hơn dùng phương pháp Runge – Kutta bậc 5: 
  i 1 i 1 3 4 5 6
16 6656 28561 9 2z y k k k k k
135 12825 56430 50 55+
= + + + − +  
Bước tính tối ưu được xác định bằng sh với s là: 
1 1
4 4
i 1 i 1 i 1 i 1
h hs 0.840896
2 z y z y+ + + +
⎛ ⎞ ⎛ ⎞ε ε= =⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟− −⎝ ⎠ ⎝ ⎠
Ta xây dựng hàm rkf() để thực hiện thuật toán trên: 
function [t, y] = rkf ( f, t0, tf, x0, parms ) 
% tim nghiem cua phuong trinh  yʹ(t) = f( t, y ), y(t0) = x0 
% dung phuong phap Runge‐Kutta‐Fehlberg bac 4, bac 5 
neqn = length(x0); 
hmin = parms(1); 
hmax = parms(2); 
tol  = parms(3); 
t(1) = t0; 
y(1:neqn, 1) = x0ʹ; 
count = 0; 
h = hmax; 
i = 2; 
while( t0 < tf ) 
    if nargin(f) > 1 
        k1 = h*feval(f, t0, x0); 
        k2 = h*feval(f, t0 + h/4, x0 + k1/4); 
        k3 = h*feval(f, t0 + 3*h/8, x0 + 3*k1/32 + 9*k2/32); 
        k4 = h*feval(f, t0 + 12*h/13, x0 + 1932*k1/2197 ‐... 
               7200*k2/2197 + 7296*k3/2197); 
        k5 = h*feval(f, t0 + h, x0 + 439*k1/216 ‐ 8*k2 +... 
               3680*k3/513 ‐ 845*k4/4104); 
        k6 = h*feval(f, t0 + h/2, x0 ‐ 8*k1/27 + 2*k2 ‐... 
               3544*k3/2565 + 1859*k4/4104 ‐ 11*k5/40); 
    else 
        k1 = h * feval(f, t0); 
        k2 = h * feval(f, t0 + h/4); 
 395
        k3 = h * feval(f, t0 + 3*h/8); 
        k4 = h * feval(f, t0 + 12*h/13); 
        k5 = h * feval(f, t0 + h); 
        k6 = h * feval(f, t0 + h/2); 
      end 
      r = max(abs(k1/360 ‐ 128*k3/4275 ‐ 2197*k4/75240 +... 
           k5/50 + 2*k6/55)/h); 
      q = 0.84*(tol/r)^(1/4); 
      if ( r < tol ) 
         x0 = x0 + 16*k1/135 + 6656*k3/12825 + 28561*k4/56430 ‐... 
                 9*k5/50 + 2*k6/55; 
          t0 = t0 + h; 
          t(i) = t0; 
          y(1:neqn, i) = x0ʹ; 
          i = i + 1;   
     end; 
      h = min(max(q, 0.1), 4.0)*h; 
     if (h > hmax) 
         h = hmax;  
      end; 
      if (t0 + h > tf) 
                    h = tf ‐ t0; 
     elseif (h < hmin) 
                   disp(ʹCan giam kich thuoc buoc tinhʹ); 
         return; 
     end; 
end; 
Để giải phương trình ta dùng chương trình ctrkf.m: 
clear all, clc 
a = 0; 
b = 1; 
y  = @f2; 
ya = 0.5; 
p = [1e‐5 1e‐3 1e‐8];%[hmin  hmax tol] 
 396
[t, y] = rkf(y, a, b, ya, p) 
plot(t, y);