CV10 - Metoda sítí - PDR-parabolické rovnice

Program ExplMet  pro výpočet řešení parciální diferenciální rovnice parabolického typu (du/dt)=(((d^2) u)/(d (u^2))) na intervalu ⟨a, b⟩ s okrajovou podmínkou u(0,t)=u(1,t)=0 a počáteční podmínkou u(x,0)=f(x) .

ExplMet[f_, m_, n_,k_] := 
  Module[{h,d1,d2,u,i,j,m1,n1}, 
	h = N[1/n]; 
	n1 = n + 1;
	m1 = m + 1; 
	alfa=N[k/(h*h)];
	Print["alfa=",alfa];
	If[alfa>0.5,Print["Metoda je nestabilni"];
	d1 = alfa; 
	  d2 = 1-2*alfa; 
	  u=Table[0,{m1},{n1}];
    	  For[i=2,i<=n,i++,u[[1,i]]=N[f[ (i-1)*h]]]; 
	  u[[1,1]] = 0.0; 
    	  u[[1,n1]] = 0.0; 
	  For[j=1,j<=m,j++,
	    For[i=2,i<=n,i++, 
	      u[[j+1,i]]=d1*(u[[j,i-1]]+u[[j,i+1]])+d2*u[[j,i]] 
            ];
	    u[[j+1,1]]=0.0;
	    u[[j+1,n1]]=0.0;
	  ],Print["Metoda je stabilni"];
	  d1 = alfa; 
	  d2 = 1-2*alfa; 
	  u=Table[0,{m1},{n1}];
    	  For[i=2,i<=n,i++,u[[1,i]]=N[f[ (i-1)*h]]]; 
	  u[[1,1]] = 0.0; 
    	  u[[1,n1]] = 0.0; 
	  For[j=1,j<=m,j++,
	    For[i=2,i<=n,i++, 
	      u[[j+1,i]]=d1*(u[[j,i-1]]+u[[j,i+1]])+d2*u[[j,i]] 
            ];
	    u[[j+1,1]]=0.0;
	    u[[j+1,n1]]=0.0;
	  ];
    	];
	res=Table[{(i-1)*h,u[[j,i]]},{j,1,m1},{i,1,n1}];
	res
   ]

Použítí programu ExplMet pro výpočet řešení parciální  diferenciální rovnice parabolického typu s počáteční podminkou f(x)=1-|1-2 x|
α=0.48

Definujeme funkci f(x)=1-|1-2 x|

f[x_] = 1 - Abs[1 - 2*x]

Použijeme program ExplMet, n=100 (dělení intervalu pro t), m=10 (dělení intervalu pro x). Krok pro t bude 0.0048. Výsledek uložíme do proměnné vys.

vys = ExplMet[f, 100, 10, 0.0048];

Zobrazíme si hodnoty  proměnné x pro t=0.0.

v0 = vys[[1]]

Zobrazíme si hodnoty  proměnné x pro t=0.048.

v10 = vys[[11]]

Zobrazíme si hodnoty  proměnné x pro t=0.096.

v20 = vys[[21]]

Zobrazíme si hodnoty  proměnné x pro t=0.192.

v40 = vys[[41]]

Zobrazíme si hodnoty  proměnné x pro t=0.48.

v100 = vys[[101]]

Grafické vykreslení řešení u(x,t) pro t=0, t=0.048, t=0.096, t= 0.192, t=0.48.

g1 = 
        ListPlot[v0, 
          PlotStyle -> {PointSize[0.02], 
              RGBColor[0, 1, 0]}];

g1l = 
        ListPlot[v0, PlotStyle -> {PointSize[0.02], RGBColor[0, 1, 0]}, 
          PlotJoined -> True];

g10 = 
        ListPlot[v10, 
          PlotStyle -> {PointSize[0.02], 
              RGBColor[1, 0, 0]}];

g10l = 
        ListPlot[v10, PlotStyle -> {PointSize[0.02], RGBColor[1, 0, 0]}, 
          PlotJoined -> True];

g20 = 
        ListPlot[v20, 
          PlotStyle -> {PointSize[0.02], 
              RGBColor[0, 0, 1]}];

g20l = 
        ListPlot[v20, PlotStyle -> {PointSize[0.02], RGBColor[0, 0, 1]}, 
          PlotJoined -> True];

g40l = 
        ListPlot[v40, PlotStyle -> {PointSize[0.02], RGBColor[1, 1, 0]}, 
          PlotJoined -> True];

g40 = 
        ListPlot[v40, 
          PlotStyle -> {PointSize[0.02], 
              RGBColor[1, 1, 0]}];

g40l = 
        ListPlot[v40, PlotStyle -> {PointSize[0.02], RGBColor[1, 1, 0]}, 
          PlotJoined -> True];

g100 = 
        ListPlot[v100, 
          PlotStyle -> {PointSize[0.02], 
              RGBColor[0, 1, 1]}];

g100l = 
        ListPlot[v100, 
          PlotStyle -> {PointSize[0.02], RGBColor[0, 1, 1]}, 
          PlotJoined -> True];

Show[{g1, g1l, g10, g10l, g20, g20l, g40, g40l, g100, 
        g100l}]

Použítí programu ExplMet pro výpočet řešení parciální  diferenciální rovnice parabolického typu s počáteční podminkou f(x)=1-|1-2 x|
α=0.52

Použijeme program ExplMet, n=100 (dělení intervalu pro t), m=10 (dělení intervalu pro x). Krok pro t bude 0.0052. Výsledek uložíme do proměnné vys.

vys = ExplMet[f, 100, 10, 0.0052];

Zobrazíme si hodnoty  proměnné x pro t=0.0.

v0 = vys[[1]]

Zobrazíme si hodnoty  proměnné x pro t=0.052.

v10 = vys[[11]]

Zobrazíme si hodnoty  proměnné x pro t=0.104.

v20 = vys[[21]]

Zobrazíme si hodnoty  proměnné x pro t=0.208.

v40 = vys[[41]]

Zobrazíme si hodnoty  proměnné x pro t=0.52.

v100 = vys[[101]]

Grafické vykreslení řešení u(x,t) pro t=0, t=0.052, t=0.104, t= 0.208, t=0.52.

g1 = 
        ListPlot[v0, 
          PlotStyle -> {PointSize[0.02], 
              RGBColor[0, 1, 0]}];

g1l = 
        ListPlot[v0, PlotStyle -> {PointSize[0.02], RGBColor[0, 1, 0]}, 
          PlotJoined -> True];

g10 = 
        ListPlot[v10, 
          PlotStyle -> {PointSize[0.02], 
              RGBColor[1, 0, 0]}];

g10l = 
        ListPlot[v10, PlotStyle -> {PointSize[0.02], RGBColor[1, 0, 0]}, 
          PlotJoined -> True];

g20 = 
        ListPlot[v20, 
          PlotStyle -> {PointSize[0.02], 
              RGBColor[0, 0, 1]}];

g20l = 
        ListPlot[v20, PlotStyle -> {PointSize[0.02], RGBColor[0, 0, 1]}, 
          PlotJoined -> True];

g40l = 
        ListPlot[v40, PlotStyle -> {PointSize[0.02], RGBColor[1, 1, 0]}, 
          PlotJoined -> True];

g40 = 
        ListPlot[v40, 
          PlotStyle -> {PointSize[0.02], 
              RGBColor[1, 1, 0]}];

g40l = 
        ListPlot[v40, PlotStyle -> {PointSize[0.02], RGBColor[1, 1, 0]}, 
          PlotJoined -> True];

g100 = 
        ListPlot[v100, 
          PlotStyle -> {PointSize[0.02], 
              RGBColor[0, 1, 1]}];

g100l = 
        ListPlot[v100, 
          PlotStyle -> {PointSize[0.02], RGBColor[0, 1, 1]}, 
          PlotJoined -> True];

Show[{g1, g1l, g10, g10l, g20, g20l, g40, g40l, g100, 
        g100l}]