CV9 - Metoda sítí - okrajové úlohy

Program MetSiti0 a MetSiti pro výpočet řešení diferenciální rovnice y''+c y'+d y = f(x) na intervalu ⟨a, b⟩ s okrajovou podmínkou y(a)=(Subscript[y, 0])  a y(b)=(Subscript[y, 1]) .

\(MetSiti0[f_, a_, b_, c_, d_, y0_, y1_, n_] := 
	
      Module[{}, 
			h = N[\((b - a)\)/n]; 
			n1 = 
          n - 1; 
			n2 = n - 2; 
			d1 = N[2 - c*h]; 
			d2 = 
          N[\(-4\) + 2*d*h*h]; 
			d3 = N[2 + c*h]; 
			A = 
          Table[0.0, {n1}, {n1}]; 
			p = 
          Table[0.0, {n1}]; 
			A[\([1, 1]\)] = 
          d2; 
			A[\([1, 2]\)] = d3; 
			A[\([n1, n2]\)] = 
          d1; 
			A[\([n1, n1]\)] = d2; 
			p[\([1]\)] = 
          2*h*h*f[a + h] - y0*d1; 
			p[\([n1]\)] = 
          2*h*h*f[a + n1*h] - y1*d3; 
			For[i = 2, i < n1, 
          A[\([i, i - 1]\)] = 
            d1; 
																			\ A[\([i, i]\)] = 
            d2; 
																			\ A[\([i, i + 1]\)] = 
            d3; 
																			\ p[\([i]\)] = 
            2*h*h*f[a + i*h]; \(i++\)]; 
			Print[
          MatrixForm[A], "   " MatrixForm[p]]; 
			p = 
          LinearSolve[A, p]; 
			res = 
          Table[{a + i*h, p[\([i]\)]}, {i, n1}]; 
			res = 
          Prepend[res, {a, y0}]; 
			res = 
          Append[res, {b, y1}]; 
			res
	]\)
\(			\)

\(MetSiti[f_, a_, b_, c_, d_, y0_, y1_, n_] := 
	
      Module[{}, 
			h = N[\((b - a)\)/n]; 
			n1 = 
          n - 1; 
			n2 = n - 2; 
			d1 = N[2 - c*h]; 
			d2 = 
          N[\(-4\) + 2*d*h*h]; 
			d3 = N[2 + c*h]; 
		\ \ A = 
          Table[Switch[i - j, \(-1\), d3, 0, d2, 1, d1, _, 0], {i, n1}, {j, 
              n1}]; 
			p = 
          Table[2*h*h*f[a + i*h], {i, n1}]; 
			p[\([1]\)] = 
          p[\([1]\)] - y0*d1; 
			p[\([n1]\)] = 
          p[\([n1]\)] - 
            y1*d3; 
 (*	\(Print[
              MatrixForm[A], "   " MatrixForm[p]];\)*) 
		\ LU = 
          LUDecomposition[A]; 
			p = 
          LUBackSubstitution[LU, p]; 
			res = 
          Table[{a + i*h, p[\([i]\)]}, {i, n1}]; 
			res = 
          Prepend[res, {a, y0}]; 
			res = 
          Append[res, {b, y1}]; 
			res		
	]\)
\(			\)

Použítí programu MetSiti0 a MetSiti pro výpočet řešení diferenciální rovnice y''+2 y'+y =(x^2) s okrajovou podmínkou y(0)=2 a y(1)=1 (c=2, d=1, počet bodů sítě n=10).

f[x_] = x*x;

data = MetSiti0[f, 0, 1, 2, 1, 6, 6, 10]

ListPlot[data, PlotJoined -> True, 
      PlotRange -> {{0, 1}, {6, 7}}]

data = MetSiti[f, 0, 1, 2, 1, 6, 6, 10]

ListPlot[data, PlotJoined -> True, 
      PlotRange -> {{0, 1}, {6, 7}}]

Použítí programu  MetSiti pro výpočet řešení diferenciální rovnice y''+2 y'+y =(x^2) s okrajovou podmínkou y(0)=2 a y(1)=1 (c=2, d=1, počet bodů sítě n=3 a n= 100, porovnání výsledků).

data1 = MetSiti[f, 0, 1, 2, 1, 6, 6, 100];
data2 = MetSiti[f, 0, 1, 2, 1, 6, 6, 3];

g1 = 
      ListPlot[data1, PlotJoined -> False, 
        PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.01]}, 
        PlotRange -> {{0, 1}, {5, 7}}]

g2 = 
      ListPlot[data2, PlotJoined -> False, 
        PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.02]}, 
        PlotRange -> {{0, 1}, {5, 7}}]

Show[{g1, g2}]

g3 = 
      ListPlot[data1, PlotJoined -> True, 
        PlotStyle -> {RGBColor[1, 0, 0]}, 
        PlotRange -> {{0, 1}, {5, 7}}]

g[x_] = InterpolatingPolynomial[data2, x]

g4 = 
      Plot[g[x], {x, 0, 1}, PlotStyle -> {RGBColor[0, 0, 1]}]

Show[{g3, g4}]

Program MetSiti3 využívá řešení soustavy s třídiagonální maticí pouze s použitím tří diagonál (neukládáme celou matici).

<<LinearAlgebra`Tridiagonal`
\(MetSiti3[f_, a_, b_, c_, d_, 
        y0_, y1_, n_] := 
	
      Module[{}, 
			h = N[\((b - a)\)/n]; 
			n1 = 
          n - 1; 
			n2 = n - 2; 
			d1 = N[2 - c*h]; 
			d2 = 
          N[\(-4\) + 2*d*h*h]; 
			d3 = N[2 + c*h]; 
			p = 
          Table[2*h*h*f[a + i*h], {i, n1}]; 
			p[\([1]\)] = 
          p[\([1]\)] - y0*d1; 
			p[\([n1]\)] = 
          p[\([n1]\)] - y1*d3; 
		\ \ \ \ diag1 = 
          Table[d1, {i, n2}]; 
		\ \ \ \ \ diag2 = 
          Table[d2, {i, n1}]; 
		\ \ \ \ \ diag3 = 
          Table[d3, {i, n2}]; 
			\ p = 
          TridiagonalSolve[diag1, diag2, diag3, p]; 
			res = 
          Table[{a + i*h, p[\([i]\)]}, {i, n1}]; 
			res = 
          Prepend[res, {a, y0}]; 
			res = 
          Append[res, {b, y1}]; 
			res		
	]\)

data = MetSiti3[f, 0, 1, 2, 1, 6, 6, 10]

g = 
      ListPlot[data, PlotJoined -> True, 
        PlotStyle -> {RGBColor[1, 0, 0]}, 
        PlotRange -> {{0, 1}, {5, 7}}]