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
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