2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




 
 Эллиптическая задача, сетка,контурный график, Mathematica
Сообщение26.06.2012, 12:07 
Здравствуйте.
Товарищи, помогите пожалуйста разобраться , что показывает выводимый 3х мерный график, после решения эллиптической задачи?
$$$-\Delta u = 0, \,\, 0<x<1, \,\, 0<y<1 $$$
$$$u(x,0)=0,   \,\,  u(x,1)=x     \,\, 0\leq x \leq 1$$$
$$$u(0,y)=0,  \,\,  u(1,y)=y,    \, \,0\leq u\leq 1$$$
Код решения задачи:

Код:
Dirichlet[n_, m_] :=
  Module[{Err = 1.0, i, j, k = 0, \[Omega]},
   \[Omega] = 4/(
    2 + Sqrt[4 - (Cos[\[Pi]/(n - 1)] + Cos[\[Pi]/(m - 1)])^2]); 
   While[ And[Err > 0.001 , k <= 70 ],
    Err = 0.0;
    For[ j = 2, j <= m - 1, j++,
     For[ i = 2, i <= n - 1, i++,
       Relax =
        1/4 (\[Omega] (u[[i, j + 1]] + u[[i, j - 1]] + u[[i + 1, j]] +
              u[[i - 1, j]] - 4 u[[i, j]]));
       u[[i, j]] = u[[i, j]] + Relax;
       Err = Max[Err, Abs[Relax]];  ]; ];
    Print["Max grid change = ", Err];
    k = k + 1;  ]; ];


n = 9;
m = 9;
Print["n = ", n]
Print["m = ", m]

u = Table[70.0, {n}, {m}];
For[i = 2, i < n, i++, u[[i, 1]] = i];
For[i = 2, i < n, i++, u[[i, m]] = 0];
For[j = 2, j < m, j++, u[[1, j]] = 0];
For[j = 2, j < m, j++, u[[n, j]] = j];

u[[1, 1]] = 1/2 (u[[1, 2]] + u[[2, 1]]);
u[[1, m]] = 1/2 (u[[1, m - 1]] + u[[2, m]]);
u[[n, 1]] = 1/2 (u[[n - 1, 1]] + u[[n, 2]]);
u[[n, m]] = 1/2 (u[[n - 1, m]] + u[[n, m - 1]]);


Dirichlet[n, m];


ListPlot3D[Transpose[u], AxesLabel -> {"y(j)", "x(i)", "u "},
ViewPoint -> {4, 2, 3}, Lighting -> False, ColorFunction -> Hue]
Print["Solução numerica para P.D.E:"]
Print["-delta u = 0"]


В Книжке Burden Analisys Numerical
http://ifolder.ru/31217077 на стр 697,701,702 есть пример решения такой задачи, но я не понял если честно также, в таблице на стр 702 (Talbe 12.2) , какой столбец соответсвутют столбцу "max greed change" из программы
http://math.fullerton.edu/mathews/n2003/ellipticpde/EllipticPDEMod/Links/EllipticPDEMod_lnk_3.html

Так же не очень понятно, зачем здесь мы задаём число 70 ? Разве гарница интервала не ограничена 1 ?
Код:
[b]u = Table[70.0, {n}, {m}];[/b]
For[i = 2, i < n, i++, u[[i, 1]] = i];
For[i = 2, i < n, i++, u[[i, m]] = 0];
For[j = 2, j < m, j++, u[[1, j]] = 0];
For[j = 2, j < m, j++, u[[n, j]] = j];

 
 
 
 Re: Эллиптическая задача, сетка,контурный график, Mathematica
Сообщение26.06.2012, 12:41 
Это приближенное решение задачи с нач. условиями $u(x,0)=x$, $u(x,1)=0$, $u(0,y)=0$, $u(1,y)=y$ :-)
Для $u(x,y)=xy$ должно быть
Код:
u = Table[70.0, {n}, {m}];
For[i = 2, i < n, i++, u[[i, 1]] = 0];
For[i = 2, i < n, i++, u[[i, m]] = i];
For[j = 2, j < m, j++, u[[1, j]] = 0];
For[j = 2, j < m, j++, u[[n, j]] = j];

А 70 здесь начальное приближение. Это типа мы решаем уравнение теплопроводности. Нач. т-ра 70, а в процессе стремится к стационарному распределению, определяемому граничными условиями.

 
 
 
 Re: Эллиптическая задача, сетка,контурный график, Mathematica
Сообщение26.06.2012, 12:47 
Спасибо :-)
Только Вы не могли бы пояснить пожалуйста, а что тут означает u=xy, точнее откуда мы это получили :-)

И как отразить в коде (и надо ли это делать) знак минус?

$$$-\Delta u = 0, \,\, 0<x<1, \,\, 0<y<1 $$$

 
 
 
 Re: Эллиптическая задача, сетка,контурный график, Mathematica
Сообщение26.06.2012, 21:06 
Поскольку уравнение однородное, то знак без разницы. А решение я просто угадал. Раз на границе полиномы, то...

 
 
 
 Re: Эллиптическая задача, сетка,контурный график, Mathematica
Сообщение27.06.2012, 00:04 
Ааа., понятно :-) Спасибо большое :-)

Вы не могли бы пожалуйста ещё пояснить, а что значат здесь значения n,m ? И почему при различных их значениях мы получаем разные графики трёхмерные ?

n=100,m=100
n=21,m=21
n=9,m=9
дают нам совершенно разные графики...

Код:
n = 100;
m = 100;
Print["n = ", n]
Print["m = ", m]

u = Table[70.0, {n}, {m}];
For[i = 2, i < n, i++, u[[i, 1]] = 0];
For[i = 2, i < n, i++, u[[i, m]] = 1];
For[j = 2, j < m, j++, u[[1, j]] = 1];
For[j = 2, j < m, j++, u[[n, j]] = 0];

u[[1, 1]] = 1/2 (u[[1, 2]] + u[[2, 1]]);
u[[1, m]] = 1/2 (u[[1, m - 1]] + u[[2, m]]);
u[[n, 1]] = 1/2 (u[[n - 1, 1]] + u[[n, 2]]);
u[[n, m]] = 1/2 (u[[n - 1, m]] + u[[n, m - 1]]);


Dirichlet[n, m];



п.с. т.е. можно сказать, что например, на графике изображено искривление металлического листа стали,задаваемое эллиптическим уравнением, верно ?

-- 27.06.2012, 01:12 --

Vince Diesel в сообщении #589238 писал(а):
Это приближенное решение задачи с нач. условиями $u(x,0)=x$, $u(x,1)=0$, $u(0,y)=0$, $u(1,y)=y$ :-)
Для $u(x,y)=xy$ должно быть
Код:
u = Table[70.0, {n}, {m}];
For[i = 2, i < n, i++, u[[i, 1]] = 0];
For[i = 2, i < n, i++, u[[i, m]] = i];
For[j = 2, j < m, j++, u[[1, j]] = 0];
For[j = 2, j < m, j++, u[[n, j]] = j];

А 70 здесь начальное приближение. Это типа мы решаем уравнение теплопроводности. Нач. т-ра 70, а в процессе стремится к стационарному распределению, определяемому граничными условиями.



Но в этом случае, получается что просто нет ограничений для x,y ?
$0\leq x \leq 1$
$0\leq y \leq 1$

-- 27.06.2012, 01:17 --

И можно у Вас ещё пожалуйста уточнить, а что изображает (охначает контурный график в этой задаче и как мы для него записываем эти параметры входные ? 0,180,15 ?
http://math.fullerton.edu/mathews/n2003/ellipticpde/EllipticPDEMod/Links/EllipticPDEMod_lnk_3.html

И как написать его же для моего номера ?

Код:
c = Table[x, {x, 0, 180, 15}];
ListContourPlot[Reverse[u], Contours -> c, ColorFunction -> Hue];
Print["Solução numerical para P.D.E."];
Print["- delta u = 0 "];

 
 
 
 Re: Эллиптическая задача, сетка,контурный график, Mathematica
Сообщение29.06.2012, 10:03 
Товарищи, ну так что Вы думаете?

 
 
 [ Сообщений: 6 ] 


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group