model.m

%Решение модельной задачи

%Параметры N1, N2 - количество узлов по осям x и t

function [U_a,U,U_n] = fun(N1,N2)

%Задание границ области дискретного изменения аргумента

L_const=1;

T_const=2^(3/2)/3;

x_begin=0;

x_end=L_const;

t_begin=0;

t_end=T_const;

%Определение шага по осям x и t

h=L_const/N1;

tau=T_const/N2;

gamma=tau/h;

%Определение координат узлов сетки

X=x_begin:h:x_end;

T=t_begin:tau:t_end;

N=length(X);

M=length(T);

%Сеточная функция внешней силы f(x,t)

for i=1:N

for j=1:M

f(i,j)=-81/4*X(i)-27/4;

end;

end;

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

%--------------Сеточный аналог теоретического решения----------------------

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

for i=1:N

for j=1:M

U_a(i,j)=27/8*X(i)^3-27/8*T(j)^2;

end;

end;

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

%-------Решение краевой задачи с применением явной разностной схемы--------

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

%Сеточная функция U(x,t)

U=zeros(N,M);

%Задание первого начального условия U(x,0)

for i=1:N

U(i,1)=27/8*X(i)^3;

end;

%Задание первого граничного условия U(L,t)

for j=1:M

U(N,j)=27/8-27/8*T(j)^2;

end

%Задание второго начального условия Ut(x,0) - по известным U(i,0)

%определяем значения U(i,1) для i=1,...,N-1

for i=2:N-1

U(i,2)=(1-gamma^2)*U(i,1)+gamma^2*U(i-1,1)/2+gamma^2*U(i+1,1)/2+tau^2*f(i,1)/2;

end

%Определение значения U(0,1) из начальных и граничных условий

U(1,2)=gamma^2*U(2,1)-(gamma^2-1)*U(1,1)+tau^2*f(1,1)/2;

%Вычисление значения функции в узлах

for j=2:M-1

%Из второго граничного условия находим U(0,j) для j=2,...,M

U(1,j+1)=2*gamma^2*U(2,j)+2*(1-gamma^2)*U(1,j)-U(1,j-1)+tau^2*f(1,j);

%Находим значения сеточной функции во внутренних узлах сетки

for i=2:N-1

U(i,j+1)=2*(1-gamma^2)*U(i,j)+gamma^2*(U(i-1,j)+U(i+1,j))-U(i,j-1)+tau^2*f(i,j);

end

end

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

%------Решение краевой задачи с применением неявной разностной схемы-------

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

%Сеточная функция U(x,t)

U_n=zeros(N,M);

%Задание первого начального условия U(x,0)

for i=1:N

U_n(i,1)=27/8*X(i)^3;

end;

%Задание первого граничного условия U(L,t)

for j=1:M

U_n(N,j)=27/8-27/8*T(j)^2;

end

%Вычисление значений сеточной функции U(i,1) методом прогонки

%Определение прогоночных коэффициентов

A=zeros(N-1,1);

B=zeros(N-1,1);

ALPHA=zeros(N-1,1);

F=zeros(N,1);

C=zeros(N,1);

BETA=zeros(N,1);

for i=1:N-2

A(i)=gamma^2/2;

B(i+1)=gamma^2/2;

C(i+1)=1+gamma^2;

F(i+1)=U_n(i+1,1)+tau^2*f(i+1,2)/2;

end;

A(N-1)=0;

B(1)=gamma^2;

C(1)=1+gamma^2;

C(N)=1;

F(1)=U_n(1,1)+tau^2*f(1,1)/2;

F(N)=27/8-27/8*T(1)^2;

%Прямой ход прогонки

ALPHA(1)=B(1)/C(1);

BETA(1)=F(1)/C(1);

for i=2:N-1

ALPHA(i)=B(i)/(C(i)-A(i-1)*ALPHA(i-1));

BETA(i)=(F(i)+A(i-1)*BETA(i-1))/(C(i)-A(i-1)*ALPHA(i-1));

end;

BETA(N)=(F(N)+A(N-1)*BETA(N-1))/(C(N)-A(N-1)*ALPHA(N-1));

%Обратный ход прогонки

U_n(N,2)=BETA(N);

for i=N-1:-1:1

U_n(i,2)=ALPHA(i)*U_n(i+1,2)+BETA(i);

end;

%Вычисление значений сеточной функции U(i,j) j=2,...,M методом прогонки

%Определение прогоночных коэффициентов

A=zeros(N-1,1);

B=zeros(N-1,1);

ALPHA=zeros(N-1,1);

F=zeros(N,1);

C=zeros(N,1);

BETA=zeros(N,1);

for i=1:N-2

A(i)=gamma^2;

B(i+1)=gamma^2;

C(i+1)=1+2*gamma^2;

end;

A(N-1)=0;

B(1)=2*gamma^2;

C(1)=1+2*gamma^2;

C(N)=1;

for j=3:M

F=zeros(N,1);

F(1)=2*U(1,j-1)-U(1,j-2)+tau^2*f(1,j-1);

F(N)=27/8-27/8*T(j)^2;

ALPHA=zeros(N-1,1);

BETA=zeros(N,1);

%Прямой ход прогонки

ALPHA(1)=B(1)/C(1);

BETA(1)=F(1)/C(1);

for i=2:N-1

F(i)=2*U(i,j-1)-U(i,j-2)+tau^2*f(i,j-1);

ALPHA(i)=B(i)/(C(i)-A(i-1)*ALPHA(i-1));

BETA(i)=(F(i)+A(i-1)*BETA(i-1))/(C(i)-A(i-1)*ALPHA(i-1));

end;

BETA(N)=(F(N)+A(N-1)*BETA(N-1))/(C(N)-A(N-1)*ALPHA(N-1));

%Обратный ход прогонки

U_n(N,j)=BETA(N);

for i=N-1:-1:1

U_n(i,j)=ALPHA(i)*U_n(i+1,j)+BETA(i);

end;

end;

Main.m

clear;

clc;

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

disp('РЕШЕНИЕ УРАВНЕНИЯ КОЛЕБАНИЙ');

disp('Решения, полученные с использованием явной и неявной схем при 50 узлах на пространственной и временной осях');

fun_graph(50,50);

[U1,U2]=fun(10,10);

disp('Таблица значений численного решения уравнения по явной схеме при 10 узлах на пространственной и временной осях');

U1

disp('Таблица значений численного решения уравнения по неявной схеме при 10 узлах на пространственной и временной осях');

U2

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

L_const=1;

disp('ИССЛЕДОВАНИЕ СХОДИМОСТИ РЕШЕНИЯ С ИСПОЛЬЗОВАНИЕМ ЯВНОЙ И НЕЯВНОЙ СХЕМ');

%Формируем вектор количества узлов по пространственной координате

Point_X=[10 20 50 100 200 500 1000];

T=1000;

%Для определения сходимости метода рассматриваем вложенные сетки

%Таблицы значений функции на самой крупной сетке

[Evident,NotEvident]=fun(Point_X(1),T);

[heigth width]=size(Evident);

last=length(Point_X);

for i=2:last

E=zeros(heigth,width);

NE=zeros(heigth,width);

mult=Point_X(i)/Point_X(1);

[Ev,NEv]=fun(Point_X(i),T);

[row column]=size(Ev);

E(1,1)=Ev(1,1);

E(1,width)=Ev(1,column);

E(heigth,1)=Ev(row,1);

E(heigth,width)=Ev(row,column);

NE(1,1)=NEv(1,1);

NE(1,width)=NEv(1,column);

NE(heigth,1)=NEv(row,1);

NE(heigth,width)=NEv(row,column);

for k=1:row

for l=1:column

if (k==1)

E(k,l)=Ev(k,l);

NE(k,l)=NEv(k,l);

end;

if (mod(k-1,mult)==0)

E((k-1)/mult+1,l)=Ev(k,l);

NE((k-1)/mult+1,l)=NEv(k,l);

end;

end;

end;

disp(strcat('Максимальное расхождение значений сеточной функции при увеличении количества узлов от',32,num2str(Point_X(i-1)),32,'до',32,num2str(Point_X(i))));

disp('Явная схема');

delta_ev=max(max(abs(Evident-E)))

disp('Неявная схема');

delta_nev=max(max(abs(NotEvident-NE)))

disp(strcat('Относительное расхождение значений сеточной функции при увеличении количества узлов от',32,num2str(Point_X(i-1)),32,'до',32,num2str(Point_X(i))));

disp('Явная схема');

%Строка, в которой находится элемент максимального расхождения

[temp index_i_ev]=max(abs(Evident-E));

%Столбец, в котором находится элемент максимального расхождения

[temp index_j_ev]=max(max(abs(Evident-E)));

d_ev=delta_ev/E(index_i_ev(index_j_ev),index_j_ev)

disp('Неявная схема');

%Строка, в которой находится элемент максимального расхождения

[temp index_i_nev]=max(abs(NotEvident-NE));

%Столбец, в котором находится элемент максимального расхождения

[temp index_j_nev]=max(max(abs(NotEvident-NE)));

d_nev=delta_nev/NE(index_i_nev(index_j_nev),index_j_nev)

Evident=E;

NotEvident=NE;

end;

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

disp('РЕШЕНИЕ МОДЕЛЬНОЙ ЗАДАЧИ');

disp('Решения, полученные с использованием явной и неявной схем при 50 узлах на пространственной и временной осях');

model_graph(50,50);

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

disp('ПРОВЕРКА СКОРОСТИ СХОДИМОСТИ РАЗНОСТНЫХ СХЕМ');

Point_X=[10 100 1000];

Point_T=[10 100 1000];

last=length(Point_X);

for i=1:last

[M_a,M,M_n]=model(Point_X(i),Point_X(i));

disp(strcat('Норма разности теоретического решения и решения по явной схеме при шаге h=',32,num2str(L_const/Point_X(i))));

d=max(max(abs(M_a-M)))

disp(strcat('Норма разности теоретического решения и решения по неявной схеме при шаге h=',32,num2str(L_const/Point_X(i))));

d_n=max(max(abs(M_a-M_n)))

end;

Список использованной литературы.

1. А.А.Самарский. Теория разностных схем. М. Наука 1983

2. В.Г. Долголаптев, В.Н. Земсков. Численные методы решения разностных уравнений математической физики. Методические указания к курсовой работе по высшей математике.- М., МИЭТ, 1987.

3. Земсков В.Н., Хахалин С.Я. Метод Сеток. Методические указания к выполнению курсовой работы на персональном компьютере. – М., МИЭТ, 1998.