Работа №3

Решить краевую задачу

.

методом конечных элементов.

Сравнить полученное решение с точным

.

Решение.

Данная задача является частным случаем вариационной задачи

,

в которой

Разделим отрезок на сегментов (конечных элементов), одинаковой длинны :

,

Решение задачи будем искать в виде

,

где

и

.

Система уравнений МКЭ будет иметь вид

,

,

где

,

,

,

,

.

Решение для

%МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ

clear all

clc

%--------------------------------------------------

a = 0;

b = 1;

%ВВОДИМ УСЛОВИЯ НА ГРАНИЦЕ

u_a = 0;

u_b = 0;

% ВВОДИМ ЧИСЛО ИНТЕРВАЛОВ

n = 4;

h = (b-a)/n; % ВЫЧИСЛЯЕМ ШАГ СЕТКИ

% ВЫЧИСЛЯЕМ КООРДИНАТЫ УЗЛОВ

X(1) = a;

for i = 1: n

X(i+1) = X(1) + i*h;

End

%СОЗДАЕМ МАТРИЦУ СИСТЕМЫ МКР И ВЕКТОР ПРАВОЙ ЧАСТИ И ЗАПОЛНЯЕМ НУЛЯМИ

A = zeros (n+1, n+1);

B = zeros (n+1, 1);

A(1, 1) = 1;

A(end, end) = 1;

B(1) = u_a;

B(end) = u_b;

for i = 2: n

r = X(i-1);

s = X(i);

t = X(i+1);

F1 = inline('-1/h^2 + (s-x).*(x-r)/h^2', 'x', 'i', 'h', 'r', 's');

A(i, i-1) = quad(F1, r, s, 1.0e-10, 0, i, h, r, s);

F2_1 = inline('1/h^2 + (x-r).^2/h^2', 'x', 'i', 'h', 'r');

F2_2 = inline('1/h^2 + (t-x).^2/h^2', 'x', 'i', 'h', 't');

A(i, i) = quad(F2_1, r, s, 1.0e-10, 0, i, h, r)+...

quad(F2_2, s, t, 1.0e-10, 0, i, h, t);

F3 = inline('-1/h^2 + (t-x).*(x-s)/h^2', 'x', 'i', 'h', 's', 't');

A(i, i+1) = quad(F3, X(i), X(i+1), 1.0e-10, 0, i, h, s, t);

F4_1 = inline('-x.*(x-r)/h', 'x', 'i', 'h', 'r');

F4_2 = inline('-x.*(t-x)/h', 'x', 'i', 'h', 't');

B(i) = quad(F4_1, r, s, 1.0e-10, 0, i, h, r)+...

quad(F4_2, s, t, 1.0e-10, 0, i, h, t);

end

%РЕШАЕМ СИСТЕМУ МКР

U = A\B

%СТРОИМ ПРИБЛИЖЕННОЕ РЕШЕНИЕ

plot(X, U, 'm-', 'LineWidth', 2)

%ВЫЧИСЛЯЕМ ТОЧНОЕ РЕШЕНИЕ

X1 = a: 0.01: b;

Y1 = exp(1)*(exp(X1)-exp(-X1))/(exp(2)-1) - X1;

%ДОБАВЛЯЕМ ГРАФИК ТОЧНОГО РЕШЕНИЯ

hold on

plot(X1, Y1, 'c-', 'LineWidth', 2)

%ОФОРМЛЯЕМ ГРАФИЧЕСКОЕ ОКНО

grid on

xlabel('{\itx}')

ylabel('{\itu(x)}')

title('-{\itu}\prime\prime + {\itu} = -{\itx}, {\itu}(0)=0 {\itu}(1)=0');

legend('МКЕ', 'ТОЧНОЕ РЕШЕНИЕ', 0)

legend('mke', 'e(e^{x}-e^{-x})/(e^2-1) - x', 0)