Метод Конечных Разностей Mathematica

[Alexeybk5]

Здравствуйте.
Товарищи, помогите пожалуйста разобраться, как можно реализовать вот такую задачу в системе Mathematica

Имеем вот такое равенство

$$u'(x_i)=\frac{1}{2h} (u(x_{i+1})-u(x_{i-1})) -\frac{h^2}{6} u'''(\xi ),$$

$ x_{i-1}<\xi <x_{i+1}, i=1,\dots,N$

Надо используя метод конечных разностей, найти систему уравнений решения задачи
$-u''(x)+p(x) u'(x)+c(x) u(x) = f(x), для x \in [a,b]$
с граничными условиями
$u(a)=\alpha$
$u(b)=\beta$

Основная сложность заключается в том, что я слабо представляю, как можно здесь работать с
$u'''(\xi ), \,\,\, x_{i-1}<\xi <x_{i+1}$
и как можно построить систему, из заданного нам равенства...


Есть решение обратной к этой задачи :
http://math.fullerton.edu/mathews/n2003/FiniteDifferenceMod.html
http://math.fullerton.edu/mathews/n2003/finitedifference/FiniteDifferenceMod/Links/FiniteDifferenceMod_lnk_2.html

$$-u''-(1/x)u'+(1+1/x)u=-(2/x)e^{x} $$

$1<x<e$
$u(1)=0 $
$u(e)=e^{e}$

Её решение выглядит таким образом :

(Оффтоп)

CoeffMat[a0_, b0_, \[Alpha]0_, \[Beta]0_, n_] :=
 Module[{a = a0, b = b0, \[Alpha] = \[Alpha]0, \[Beta] = \[Beta]0, h, 
   j},
  h = (b - a)/n; 
  Vx = Table[ a + h j, {j, n - 1} ]; 
  Va = Vc = Table[ 0, {n - 2} ]; 
  Vd = Vb = Table[ 0, {n - 1} ]; 
  
  For[ j = 1, j <= n - 1, j++,
   Vd[[j]] = 2 + h^2 q[Vx[[j]]]; 
   Vb[[j]] = - h^2 r[Vx[[j]]];  ]; 
  
  Vb[[1]] = Vb[[1]] + (1 + 1/2 h p[Vx[[1]]]) \[Alpha]; 
  Vb[[n - 1]] = Vb[[n - 1]] + (1 - 1/2 h p[Vx[[n - 1]]]) \[Beta]; 
  
  
  For[ j = 1, j <= n - 2, j++,
   Va[[j]] = -1 - 1/2 h p[Vx[[j + 1]]];
   Vc[[j]] = -1 + 1/2 h p[Vx[[j]]];  ];  ]  
(************************************************************)


TriMat[va_, vd_, vc_, vb_, n_] :=
 Module[{a = va, b = vb, c = vc, d = vd, k, u}, 
  u = Table[0, {n - 1}]; 
  
  For[ k = 2, k <= n - 1, k++, 
   d[[k]] = d[[k]] -  a[[k - 1]]/d[[k - 1]]  c[[k - 1]]; 
   b[[k]] = b[[k]] -  a[[k - 1]]/d[[k - 1]]  b[[k - 1]]; ]; 
  u[[n - 1]] = b[[n - 1]]/d[[n - 1]]; 
  
  For[ k = n - 2, 1 <= k, k--,
   u[[k]] = (b[[k]] - c[[k]] u[[k + 1]])/d[[k]] ]; 
  Return[u]; ]  
(*************************************************************)

FineDDE[a0_, b0_, \[Alpha]0_, \[Beta]0_, n_] :=
 Module[{a = a0, b = b0, \[Alpha] = \[Alpha]0, \[Beta] = \[Beta]0, i},
  CoeffMat[a, b, \[Alpha], \[Beta], n]; 
  Xt = TriMat[Va, Vd, Vc, Vb, n]; 
  pts = Transpose[{Vx, Xt}];  
  pts = Append[pts, {b, \[Beta]}];  
  pts = Prepend[pts, {a, \[Alpha]}];  
  Return[pts]; ]  
(****************************************************************)
(****************************************************************)
(****************************************************************)


p[x_] = -1/x;
q[x_] = 1 + 1/x;
r[x_] = (2/x)*Exp[x];
a = 1.0;
b = Exp[1];
\[Alpha] = 0;
\[Beta] = Exp[Exp[1]];

Print["Vamos solucar o D.E.:"]
Print["-u'' + (", p[x], ")u' +(", q[x], ")u = ", -r[x]];
Print["com u(", a, ") = \[Alpha] = ", \[Alpha]];
Print["e u(", b, ") = \[Beta] = ", \[Beta]];
Print["e intevalo x em [", a, ",", b, "]"];

(*********************************************************************)

CoeffMat[1.0, Exp[1], 0.0, Exp[Exp[1]], 50]




Print["\n", "Vx = ", Vx];
Print["\n", "Va = ", Va];
Print["\n", "Vc = ", Vc];
Print["\n", "Vd = ", Vd];
Print["\n", "Vb = ", Vb];

Xt = TriMat[Va, Vd, Vc, Vb, 50];
Print["\n", "Vx = ", Vx];
Print["\n", "Xt = ", Xt];

wi = Xt;
wi = Append[wi, 0];
wi = Prepend[wi, 1];

pts = Transpose[{Vx, Xt}];
Print["\n", "pontos = ", pts];

pts = Append[pts, {1, 0}];
Print["\n", "pontos com boundory direita: ", pts];

pts = Prepend[pts, {0, 1}];
Print["\n", "pontos com boundorry esqueda : ", pts];

pontos4 = FineDDE[1.0, Exp[1], 0.0, Exp[Exp[1]], 50];
 
Print[""]
Print["Grafico de pontos calculados de metodo das diferenças finitas:"]
Print[ListLinePlot[pontos4, PlotRange -> {{-5, 5.1}, {-5, 5.1}}, 
  AxesLabel -> {"x", "u"}]]
Print["grafico de função exacta -u'' -(1/x) * u = -(2/x) e^{x}:"]
(*Print[Plot[(1-x)*Exp[x],{x,0,1}]]*)
Print["u(", a, ") = ", \[Alpha]];
Print["u(", b, ") = ", \[Beta]];

[TOTAL]

Alexeybk5 в сообщении #589246 писал(а):

Надо используя метод конечных разностей, найти систему уравнений решения задачи
$-u''(x)+p(x) u'(x)+c(x) u(x) = f(x), для x \in [a,b]$
с граничными условиями
$u(a)=\alpha$
$u(b)=\beta$

Аппроксимируйте производные конечными разностями.

[Alexeybk5]

Да, но что делать с производной третьего порядка ? Она же тут в качестве погрешности выступает, разве нет ?

-- 26.06.2012, 14:30 --

И на каких точках можно тогда производить апроксимацию ? Т.е. может быть здесь неполное условие ?

[ewert]

Alexeybk5 в сообщении #589252 писал(а):

Да, но что делать с производной третьего порядка ? Она же тут в качестве погрешности выступает, разве нет ?

Да, поэтому её надо просто выкинуть (этот член лишь указывает на порядок точности схемы).

[Alexeybk5]

ааа... Т.е. получается, надо просто написать программу, которая будет для любых входных границах$ [a,b],\,\,\,\,\, u(a) = \alpha, \,\,\,\,\,\, u(b) = \betta$
Разбивать этот отрезок, положим на 50 частей (n=50)
и по методу конечных разностей найти
$u'(x)$
$u''(x)$

?


Но в программе так же должно быть заложена возможность менять значения альфа и бетта ? Т
$u(a) = \alpha, \,\,\,\,\,\, u(b) = \betta$


т.е. в итоге, надо написать функцию с входными параметрами такова вида:
$DiffMethod[a,b,\alpha,\betta,n]$
и программа должна выводить
$u'(x) = ...$
$u''(x) =...$

где корень этой функции будет формула подсчёта u'', u' (http://alexandr4784.narod.ru/B13/b13_7_72.pdf)
Но я всё равно не понимаю, как тогда находить эти самые $u_i$ ? задавая только $x_i$ ?

[Alexeybk5]

В работе эта задача идёт третьей, но в первом номере работы задана функция $u(x)=(1-x)e^x$ ( но совершенно не понятно, можно ли её использовать?

Но я если честно не очень понимаю само условие задачи:


"установить (найти) систему уравнений, соответствуюзую решению задачи":

$$-u''(x)+p(x)u'(x)+c(x)u(x) = f(x)$$

для$ x \in [a,b]$ и с граничными условиями
$u(a) = \alpha$
$u(b)=\beta$

[Alexeybk5]

Здравствуйте, будьте добры, если не сложно, удалите пожалуйста эту тему, я не правильно задал вопрос, сейчас перечитав задание, я во всём разобрался, и не хотелось бы засорять форум не нужными темами, спасибо