ROBNI PROBLEMI

Uvod v numeriène metode

Bor Plestenjak

8.1.1995

Diferenèna metoda

Uvod

Pri diferenèni metodi interval [a,b] razdelimo z delilnimi toèkami 

		x_0=a, x_1=a+h, ... , x_n=b, 	h=(b-a)/h.

Odvode nadomestimo s simetriènimi diferencami:

		        y_(j+1)-2 y_j+y_(j-1)				y_(j+1)-y_(j-1)
	y''(x_j) = -------------------------,          y'(x_j) = --------------------,          j=1,...,n-1
				h^2							2h

Tako dobimo n-1 enaèb 

		         y_(j+1)-2 y_j+y_(j-1)				y_(j+1)-y_(j-1)
	(*)	       ------------------------ =  f( x_j,  y_j,   -------------------- ),       j=1,...,n-1, 
				h^2							2h		

iz robnih pogojev pa dobimo še dve enaèbi:   y_0=A  in  y_n=B.

Èe je robni problem linearen, je dobljeni sistem tridiagonalen linearni sistem in ga rešimo z eno od metod za reševanje tridiagonalnih sistemov. Drugaèe pa imamo veè dela z reševanjem nelinearnega sistema. V našem primeru bomo diferenèno metodo uporabili le na linearnih robnih problemih.

Robni pogoji

V robnih pogojih lahko nastopajo tudi odvodi. Splošna robna pogoja imata obliko:

	Ay(a) + Cy'(a) = E
	By(b) + Dy'(b) = F.

Problem rešimo tako, da vpeljemo novi neznanki y_(-1) in y_(n+1), nato pa tudi odvode v robnih pogojih zapišemo s simetriènimi diferencami. Enaèbe (*) moramo sedaj zapisati v toèkah j=0,..n, saj vrednosti na robu niso znane. Tako dobimo n+3 enaèb za n+3 neznank, sistem pa je spet tridiagonalen.

Program

Program za reševanje robnega problema  iz dobljene enaèbe in robnih pogojev sestavi tridiagonalen sistem, ga reši in vrne vrednosti y_0 do y_n. Seveda je mono spremeniti število delilnih toèk n in na ta naèin spreminjati natanènost raèunanja.

Robni::usage="Robni[en,{ra,rb},a,b,n] vrne tabelo parov {x,y(x)} 
numeriène\n rešitve linearnega robnega problema, podanega z enaèbo en\n in 
robnima problemoma ra in rb. Rešitev se išèe na intervalu [a,b]\n po 
diferenèni metodi z razdelitvijo intervala na n toèk.\n Pravilna uporaba je 
naprimer:\n\n Robni[y''[x]+y[x]-x^2==0,{y[1]==0,y[2]-y'[2]==0},1,2,10]";

Robni[en_,ro_,a_,b_,n_]:=Module[

    {(*s,h,g,roba,robb,spmeja,zgmeja,enacbe,sistem,neznanke,resitev,rez*)},

	h:=N[(b-a)/n];
	g[z_]:=en /. x->z;
	roba=ro[[1]] /. y[j_]->y[z[0]];
	roba=roba /. y'[j_]->(y[z[1]]-y[z[-1]])/(2h);
	spmeja=If[StringPosition[ToString[roba],"y[z[-1]]"]=={},1,0];
	robb=ro[[2]] /. y[j_]->y[z[n]];
	robb=robb /. y'[j_]->(y[z[n+1]]-y[z[n-1]])/(2h);
	zgmeja=If[StringPosition[ToString[robb],ToString[y[z[n+1]]]]=={},n-1,n];
	enacbe=Table[g[z[i]],{i,spmeja,zgmeja}];
	enacbe=enacbe /. y'[z[j_]]->(y[z[j+1]]-y[z[j-1]])/(2h);
	enacbe=enacbe /. y''[z[j_]]->(y[z[j+1]]-2y[z[j]]+y[z[j-1]])/h^2;
	sistem=Union[enacbe,{roba,robb}];
	sistem=sistem /. y[z[j_]] -> s[j];
	sistem=sistem /. z[j_] -> a+ j h;
	neznanke=Table[s[j],{j,spmeja-1,zgmeja+1}];
	resitev=First[NSolve[sistem,neznanke]];
	rez=Table[{N[a+j h],s[j]},{j,spmeja-1,zgmeja+1}] /. resitev;
	If[spmeja==0,rez=Drop[rez,1]];
	If[zgmeja==n,rez=Drop[rez,-1]]; 
	rez ];

Zgled

Zgled je robni problem:   y''-3y'+2y=2,   y(0)=1,   y(1)=4

Toèna rešitev tega problema je 

					     2x	       x
					3(e     -   e   )
			y(x) = 1 + ----------------
					      2	
					    e	 -   e

Lahko narišemo tudi graf toène rešitve:

tocna1=Plot[1+3(E^(2x)-E^x)/(E^2-E),{x,0,1}]

Za numerièno rešitev uporabimo ukaz Robni. Navodila za uporabo so:

?Robni

pari1=Robni[y''[x]-3y'[x]+2y[x]==2,{y[0]==1,y[1]==4},0,1,4]

Kot rezultat dobimo pare toèk (x,y(x)), iz katerih lahko narišemo graf:

pribl1=ListPlot[pari1,PlotJoined->True]

Dobljeni graf lahko primerjamo s toèno rešitvijo:

Show[tocna1,pribl1]

S poveèanjem števila toèk dobimo bolj toène rezultate, a poveèamo raèunski èas.

pari2=Robni[y''[x]-3y'[x]+2y[x]==2,{y[0]==1,y[1]==4},0,1,20]

pribl2=ListPlot[pari2,PlotJoined->True]

Show[tocna1,pribl2]

Strelska metoda

Uvod

Pri strelski metodi robni problem 

(1)	y'' = f ( x, y, y' ),     y(a)=A,    y(b)=B,	

rešujemo tako, da rešujemo zaèetni problem

(2)	y'' = f ( x, y, y' ),     y(a)=A,    y'(a)=C.	

Parameter C je neznan in ga moramo doloèiti tako, da bo rešitev usrezala pogoju y(b)=B.

Ko z eno od numeriènih metod za reševanje zaèetnega problema rešimo problem (2), dobimo rešitev y(x,C), saj je rešitev odvisna od C. Rešiti moramo enaèbo y(b,C)=B.

Zgledi

Zgled 1

Za prvi zgled vzemimo enaèbo, katere rešitev je Legendrov polinom reda 4:

(x^2-1) y''+ 2x y'-12y=0,

z robnimi pogoji

y(0)=0
y(1)=1

Toèna rešitev je: y(x)=1/2(5x^3-3x), kar pomeni y'(0)=-3/2.

Pri reševanju s strelsko metodo rešujemo enaèbo:

(x^2-1) y''+ 2x y'-12y=0,

z zaèetnimi pogoji

y(0)=0
y'(0)=a

Poglejmo si, kaj dobimo pri parametrih a= -1.2 in a= -1.8.

a= - 1.2

res1=NDSolve[
	{(x^2-1) y''[x]+2 x y'[x]-12 y[x]==0,
	  y[0]==0,y'[0]==-1.2}, y[x],{x,0,1.2}]

Plot[y[x] /. res1,{x,0,1},PlotRange->All]

res1 /. x->1

a= - 1.8

res2=NDSolve[
	{(x^2-1) y''[x]+2 x y'[x]-12 y[x]==0,
	  y[0]==0,y'[0]==-1.8}, y[x],{x,0,1.2}]

Plot[y[x] /. res2,{x,0,1},PlotRange->All]

res2 /. x->1

Pravo rešitev bi dobili, ko bi izbrali a= -3/2.

res3=NDSolve[
	{(x^2-1) y''[x]+2 x y'[x]-12 y[x]==0,
	  y[0]==0,y'[0]==-1.5}, y[x],{x,0,1.1},AccuracyGoal->10]

res3 /. x->1

Zgled 2

Za drugi zgled vzemimo streljanje s topom. Na izstreljeno zrno med letom delujeta sila tee in sila zraènega upora. Poenostavljeni enaèbi sta:

x''(t) = - 0.1 Sqrt[x'[t]^2+y'[t]^2] x'[t],
y''(t) = - 0.1 Sqrt[x'[t]^2+y'[t]^2] y'[t] - 9.8

Robni pogoji naj bodo:

x(0) = 0	x(10)=200
y(0) = 0	y(10)=0

Vprašanje je, s kakšno zaèetno hitrostjo in pod kakšnim kotom je treba izstreliti zrno, da bo to izpolnjeno.

Zrno izstrelimo z zaèetno hitrostjo v in kotom a:

x'(0) = v cos(a)
y'(0) = v sin(a)

Resitev[v_,a_]:=Module[{},
	nds=NDSolve[{x''[t]==-0.01 Sqrt[x'[t]^2+y'[t]^2] x'[t],
		     y''[t]==-0.01 Sqrt[x'[t]^2+y'[t]^2] y'[t]-9.8,
		 	 x[0]==0,y[0]==0,
		     x'[0]==v Cos[a],y'[0]==v Sin[a]},
			{x,y},{t,0,10}];
	{p,q}={x,y} /. First[nds];
	Print["x[10]=",p[10]];
	Print["y[10]=",q[10]];
	ParametricPlot[{p[t],q[t]},{t,0,10},PlotRange->All] ]

Resitev[200,0.75]

Resitev[180,1]

Resitev[180,0.9]

Resitev[190,0.885]