Вот здесь уже есть немало примеров на реализацию численных методов в Excel.

Однако, решать в Excel краевую задачу для дифференциального уравнения второго порядка, по-моему, извращение. Гораздо проще написать программку на Паскале, чем мы сейчас и займёмся. Но сначала немного теории, только самые нужные формулы.

Для дифференциального уравнения второго порядка

поставлена краевая задача:
Здесь u(t) -искомое решение, A, B - краевые условия,заданные функции коэффициентов. Введем на отрезке [a,b] разностную сетку (t₀, t₁, t₂, ..., t_M,), t_i=a+τ*i, i=0,1,...,M, τ=(b-a)/M, M – параметр задачи (число шагов сетки). Заменим уравнение разностной схемой

Собирая коэффициенты при y_i, y_i+1 , y_i-1, получим следующую СЛАУ на вектор неизвестных (y₀, y₁, y₂, ..., y_M):

(k₁—k₂/τ)y₀ + k₂ y₁/τ =A

AA_i y_i-1 - CC_iy_i + BB_iy_i+1 = FF_i, i=1,2,...,M-1

(l₁+l₂/τ) y_M - l₂ y_M-1 /τ=B,

где CC_i = 2 - g(t_i) τ² , AA_i = 1 - p(t_i) τ/2, BB_i = 1 + p(t_i) τ/2, FF_i = τ²f(t_i).

Система решается методом прогонки, в котором решение отыскивается в виде y_i=α_i+1y_i+1+β_i+1, i=0,1,…,M-1. Алгоритм состоит из следующих шагов.

1. Поскольку y₀=α₁y₁+β₁, и из первого уравнения y₀ = (A - k₂ y₁/τ) / (k₁—k₂/τ), то α₁ = - k₂/ [ τ (k₁—k₂/τ) ], β₁ = A/(k₁—k₂/τ)

2. Следующие прогоночные коэффициенты a_i, b_i определяются по рекурентным формулам (прямой ход метода прогонки с трёхдиагональной матрицей):
α_i+1=, β _i+1= i=1,2,...,M-1.

3. Для определения решения на правой границе отрезка воспользуемся последним уравнением системы (l₁+l₂/τ) y_M - l₂y_M-1 /τ=B и соотношением y_M-1=α_M y_M + β_M, откуда y_M определяется как решение системы двух линейных алгебраических уравнений.

4. Обратный ход метода прогонки для определения решения.

y_i=α_i+1y_i+1+β_i+1, i = M-1,M-2, ..., 0.

Теперь напишем программу решения краевой задачи для ОДУ второго порядка с помощью конечно-разностного метода. Это программа будет тестовой, то есть, содержит функцию yt(t), определяющую точное решение.

Уравнение:

Точное решение:
u_t(t) = cos(sin(t)) + C₁ sin(sin(t)), C₁=(10 - cos(sin(1)))/sin(sin(1)).

Программа для заданного константой M будет выдавать на экран таблицу со столбцами
i t_i y_i u(t_i) |u(t_i)-y_i|,
также по последнему столбцу определяется максимальная погрешность.

Метод показывает первый порядок аппроксимации, то есть, при увеличении числа узлов в 2 раза точность решения также повышается примерно в 2 раза. Возможно, здесь чуть-чуть иначе вычисляются α1 и β1, но с посчитанной ниже задачей в Excel совпало.

program odu2test;
{Решение тестовой задачи с известным точным решением yt(t)}
const n=10; {Число точек}
      a=0;  {Границы интервала}
      b=1;
      {Краевые условия}
      k1=1; k2=0; a1=1;
      l1=1; l2=0; b1=10;
type vec=array [0..n] of real;
var
 h,eps,epsmax:real;
 al,bet,t,u,aa,bb,cc,ff:vec;
 i,j:integer;

{Функция коэффициентов при 1-й производной u'(t)}
function p(t:real):real;
begin
 p:=sin(t)/cos(t);
end;
{Функция коэффициентов при u(t)}
function g(t:real):real;
begin
 g:=sqr(cos(t));
end;
{Функция правой части уравнения}
function f(t:real):real;
begin
 f:=0;
end;
{Точное решение yt(t)}
function yt(t:real):real;
var c1:real;
begin
 c1:=(10-cos(sin(1)))/sin(sin(1));
 yt:=cos(sin(t))+c1*sin(sin(t));
end;

begin
 {Вычисление шага сетки}
 h:=(b-a)/n;
 {Формирование векторов коэффициентов для 3-диагональной СЛАУ}
 for i:=0 to n do begin
  t[i]:=a+h*i;
  aa[i]:=1-p(t[i])*h/2;
  bb[i]:=1+p(t[i])*h/2;
  cc[i]:=2-g(t[i])*h*h;
  ff[i]:=h*h*f(t[i]);
 end;
 {Краевое условие при t=a}
 al[1]:=k2/(k2-k1*h);
 bet[1]:=-(a1*h)/(k2-k1*h);
 {Прямой ход метода прогонки}
 for i:=1 to n-1 do begin
  al[i+1]:=bb[i]/(CC[i]-al[i]*aa[i]);
  bet[i+1]:=(aa[i]*bet[i]-ff[i])/(cc[i]-al[i]*aa[i]);
 end;
 {Решение на правой границе отрезка}
 u[n]:=(l2*bet[n]+b1*h)/(l2+h*l1-l2*al[n]);
 {Обратный ход метода прогонки}
 for i:=n-1 downto 0 do
  u[i]:=al[i+1]*u[i+1]+bet[i+1];
 {Вывод таблицы и определение макс.погрешности}
 writeln('i':4,'t[i]':10,'y[i]':10,'yt(t[i])':10,'Погрешн.':10);
 for i:=0 to n do begin
  eps:=abs(u[i]-yt(t[i]));
  if eps>epsmax then epsmax:=eps;
  writeln (i:4,t[i]:10:4,u[i]:10:4,yt(t[i]):10:4,eps);
 end;
 writeln ('Максимальная погрешность=',epsmax);
 reset (input); readln;
end.

Используя эту программу, можно решить и другую аналогичную краевую задачу, изменив в нашей программке всего несколько строчек. Например, такую:
Уравнение , интервал [1.1,1.4], краевое условие при t=a имеет вид u-2u ′ =-1, краевое условие при t=b имеет вид u=4. Вроде бы, должно получиться так:

program odu2my;
const n=10; {Число точек}
      a=1.1;  {Границы интервала}
      b=1.4;
      {Краевые условия}
      k1=1; k2=-2; a1=-1;
      l1=1; l2=0; b1=4;
type vec=array [0..n] of real;
var
 h:real;
 al,bet,t,u,aa,bb,cc,ff:vec;
 i,j:integer;

{Функция коэффициентов при 1-й производной u'(t)}
function p(t:real):real;
begin
 p:=-1;
end;
{Функция коэффициентов при u(t)}
function g(t:real):real;
begin
 g:=2/t;
end;
{Функция правой части уравнения}
function f(t:real):real;
begin
 f:=t+0.4;
end;

begin
 {Вычисление шага сетки}
 h:=(b-a)/n;
 {Формирование векторов коэффициентов для 3-диагональной СЛАУ}
 for i:=0 to n do begin
  t[i]:=a+h*i;
  aa[i]:=1-p(t[i])*h/2;
  bb[i]:=1+p(t[i])*h/2;
  cc[i]:=2-g(t[i])*h*h;
  ff[i]:=h*h*f(t[i]);
 end;
 {Краевое условие при t=a}
 al[1]:=k2/(k2-k1*h);
 bet[1]:=-(a1*h)/(k2-k1*h);
 {Прямой ход метода прогонки}
 for i:=1 to n-1 do begin
  al[i+1]:=bb[i]/(CC[i]-al[i]*aa[i]);
  bet[i+1]:=(aa[i]*bet[i]-ff[i])/(cc[i]-al[i]*aa[i]);
 end;
 {Решение на правой границе отрезка}
 u[n]:=(l2*bet[n]+b1*h)/(l2+h*l1-l2*al[n]);
 {Обратный ход метода прогонки}
 for i:=n-1 downto 0 do
  u[i]:=al[i+1]*u[i+1]+bet[i+1];
 {Вывод таблицы и определение макс.погрешности}
 writeln('i':4,'t[i]':10,'y[i]':10);
 for i:=0 to n do writeln (i:4,t[i]:10:4,u[i]:10:4);
 reset (input); readln;
end.

Ниже то же самое посчитано в Excel для 10, 20 и 40 узлов сетки, в файле серым фоном обозначены ячейки с формулами и константами, которые являются входными данными, жёлтым - формулы, которые считаются иначе, чем соседние формулы выше и/или ниже.

 odu2.xls (56 Кб, Excel XP/2003)

гостевая; E-mail