Matematyka.pl

Metoda Strzału

Rozpatrujemy zagadnienie brzegowe:

\(\displaystyle{ \begin{cases} y''+ p(x)y' +q(x)y = r(x) \\ y(a) = \alpha, \ \ y(b)=\beta \ \ (1) \end{cases} }\)

Rozwiązanie zagadnienia \(\displaystyle{ (1) }\) Metodą Strzału polega na zamianie zagadnienia brzegowego na kombinację liniową dwóch zagadnień początkowych drugiego rzędu:

\(\displaystyle{ u''+ p(x)u' +q(x)u = r(x) \\ u(a) = \alpha, \ \ u'(\alpha)= 0 \ \ (2)}\)

\(\displaystyle{ v''+ p(x)v' +q(x)v = r(x) \\ v(a) = 0, \ \ v'(a)= 1 \ \ (3) }\)

\(\displaystyle{ y(x) = u(x)+A v(x) \ \ (4) }\)

Na podstawie \(\displaystyle{ (2), (3), (4), (1) }\)

\(\displaystyle{ u''(x)+ Av''(x)+p(x)[u'(x)+Av"(x)] + q(x)[u(x)+ Av(x)]=r(x) }\)

\(\displaystyle{ u''(x) +p(x)u'(x)+q(x)u(x) +A[v''(x)+p(x)v'(x)+q(x)v(x)]=r(x) }\)

Stałą \(\displaystyle{ A }\) znajdujemy z warunków brzegowych:

\(\displaystyle{ y(a) = u(a)+Av(a) = \alpha +0 = \alpha, }\)

\(\displaystyle{ y(b) = u(b) +Av(b) = \beta }\)

\(\displaystyle{ A = \frac{\beta - u(b)}{v(b)}.}\)

Jeśli \(\displaystyle{ v(b) \neq 0 }\) to rozwiązanie zagadnienia brzegowego \(\displaystyle{ (1) }\) ma postać:

\(\displaystyle{ y(x) = u(x) + \frac{\beta - u(b)}{v(b)}v(x).}\)

Do rozwiązań zagadnień początkowych \(\displaystyle{ (2), (3) }\) stosujemy Metodę Runge-Kutta 4 rzędu.

Z warunkami istnienia rozwiązań tej metody można zapoznać się w podręczniku:

DAVID KINCAID WARD CHENEY analiza numeryczna. WNT Warszawa 2006.

Program w MATLAB

Kod: Zaznacz cały

 
function mstrzal(f1,f2,a,b,alfa,beta,n) 
%Rozwiązanie zagadnienia początkowego drugiego rzędu Metodą Strzałów. 
% f1= -p(x)*u-q(x)*v+r(x), f2= -p(x)*u-q(x)*v. 
h=(b-a)/n; 
y1=alfa; 
y2=0; 
u=0; 
v=alfa; 
for i=1:n 
% Runge-Kutta rzędu 4 
 x=a+(i-1)*h; 
 k1=feval(f1,x,u,v); 
 c1=u; 
 k2=feval(f1,x+h/2,u+h*k1/2,v+h*c1/2); 
 c2=u+h/2*k1; 
 k3=feval(f1,x+h/2,u+h*k2/2,v+h*c2/2); 
 c3=u+h/2*k2; 
 k4=feval(f1,x+h,u+h*k3,v+h*c3); 
 c4=u+h*k3; 
 u=u+h*(k1+2*k2+2*k3+k4)/6; 
 v=v+h*(c1+2*c2+2*c3+c4)/6; 
 y1(i+1)=v; 
end 
u=1; 
v=0; 
for i=1:n 
% Runge-Kutta rzędu 4 
 x=a+(i-1)*h; 
 k1=feval(f2,x,u,v); 
 c1=u; 
 k2=feval(f2,x+h/2,u+h*k1/2,v+h*c1/2); 
 c2=u+h/2*k1; 
 k3=feval(f2,x+h/2,u+h*k2/2,v+h*c2/2); 
 c3=u+h/2*k2; 
 k4=feval(f2,x+h,u+h*k3,v+h*c3); 
 c4=u+h*k3; 
 u=u+h*(k1+2*k2+2*k3+k4)/6; 
 v=v+h*(c1+2*c2+2*c3+c4)/6; 
 y2(i+1)=v; 
 end 
fprintf('\n') 
disp(' Metoda Strzału ') 
fprintf('\n') 
disp([' u(b) = ',num2str(y1(n+1)),' ']) 
disp([' v(b) = ',num2str(y2(n+1)),' ']) 
fprintf('\n') 
disp('_______________________________________________') 
disp('   xi         ui          vi           yi      ') 
disp('_______________________________________________') 
for i=1:n+1 
 x=a+(i-1)*h; 
 w=y1(i)+(beta-y1(n+1))/y2(n+1)*y2(i); 
  s='n'; 
 if (s=='n') 
 fprintf('%6.2f %12.6f %12.6f %12.6f\n',x,y1(i),y2(i),w) 
 else 
 err=abs(w-s); 
 fprintf('%6.2f %12.6f %12.6f %12.6f %12.6f %10.2e\n', x, y1(i),y2(i), w, s, err) 
 end 
end

Przykład

Kod: Zaznacz cały

function f=f1(x,u,v)
f=3/x*u-3/x.^2*v+2*x.^2*exp(x);
 
function f=f2(x,u,v)
f=3/x*u-3/x.^2*v;

mstrzal('f1','f2',1,2,0,4*exp(2),10)

 Metoda Strzału

 u(b) = 13.2463 
 v(b) = 3 

_______________________________________________
   xi         ui          vi           yi      
_______________________________________________
  1.00     0.000000     0.000000     0.000000
  1.10     0.032983     0.115500     0.660913
  1.20     0.158380     0.264001     1.593649
  1.30     0.423715     0.448502     2.862042
  1.40     0.888400     0.672003     4.541813
  1.50     1.625684     0.937504     6.722522
  1.60     2.724900     1.248006     9.509811
  1.70     4.294042     1.606508    13.027984
  1.80     6.462735     2.016010    17.422977
  1.90     9.385635     2.479512    22.865754
  2.00    13.246340     3.000015    29.556224