Форум кафедры Техники и Электрофизики Высоких Напряжений

Онлайн-сообщество ТВНщиков
Гостям форума:

Добро пожаловать на форум по технике высоких напряжений!
Для получения доступа ко всем разделам необходимо зарегистрироваться


Текущее время: 16 дек 2018, 11:48

Часовой пояс: UTC + 3 часа




Начать новую тему Ответить на тему  [ Сообщений: 2 ] 
Автор Сообщение
СообщениеДобавлено: 22 дек 2013, 14:31 
Не в сети

Зарегистрирован: 13 сен 2011, 12:23
Сообщения: 46
Откуда: Тамбов
Здравствуйте, Даниил Анатольевич! Не могли бы вы помочь найти ошибку? Я просмотрел графики, вроде пролонгирует и интерполирует нормально.
clear, close all, clc
A = 200; %параметр распределения заряда
B = 0.1; %параметр распределения заряда
% Параметры грубой сетки
k = 5;
l = 10;
Ngrids = 5;
maxdelta = 1;
for i = 1 : Ngrids
Xlims(i) = 2^(i-1)*l;
Ylims(i)= 2^(i-1)*k;
end


%Начальные условия
Fi = zeros(Ylims(end)+1,Xlims(end)+1);
%Граничные условия
Fi(1,:) = 0;
Fi(end,:) = 0;
Fi(:,1) = 0;
Fi(:,end) = 0;

% Параметр выхода из цикла
while maxdelta > 10^(-3)

% Расчет f в точках X,Y
x = linspace(0,2,Xlims(end)+1);
y = linspace(0,1,Ylims(end)+1);
delta = x(2)-x(1);

[X Y] = meshgrid(x,y);

f = -B*exp(-A*((X-x(end)/2).^2+(Y-y(end)/2).^2));

%Делаем одну итерацию метода Якоби
[M N] = size(Fi);

for m = 2:M-1
for n = 2:N-1
Fi(m,n) = 1/4*(Fi(m-1,n)+ Fi(m+1,n) + Fi(m,n - 1) + Fi(m,n+1)) - delta^2 /4 * f(m,n);
end
end

D{1} = zeros(M,N);
%Находим невязку
for m = 2:M-1
for n = 2:N-1
D{1}(m,n) = (Fi(m-1,n)-2*Fi(m,n)+Fi(m+1,n))/delta^2 + (Fi(m,n-1)-2*Fi(m,n)+Fi(m,n+1))/delta^2 - f(m,n);
end
end

for i = 2: Ngrids
% Задаем матрицу для невязки на i-м сеточном уровне
D{i} = zeros(Ylims(end-i+1)+1, Xlims(end-i+1)+1);

[M N] = size(D{i});
%Пролонгируем невязку на более грубую сетку
for m = 2:M-1
for n = 2:N-1
D{i}(m,n) = D{i-1}(m*2-1,n*2-1);
end
end

% Расчет f в точках X,Y
x = linspace(0,2,N);
y = linspace(0,1,M);
delta = x(2)-x(1);
[X Y] = meshgrid(x,y);
f = -B*exp(-A*((X-x(end)/2).^2+(Y-y(end)/2).^2));
% Прогоняем невязку через метод Якоби (1 итерация)
for m = 2:M-1
for n = 2:N-1
D{i}(m,n) = 1/4*(D{i}(m-1,n)+ D{i}(m+1,n) + D{i}(m,n - 1) + D{i}(m,n+1)) - delta^2 /4 * f(m,n);
end
end

end

for i = 1: Ngrids-1
% Интерполируем невязку на более мягкую сетку
x = linspace(0,2,Xlims(i)+1);
y = linspace(0,1,Ylims(i)+1);
[X Y] = meshgrid(x,y);
xx = linspace(0,2,Xlims(i+1)+1);
yy = linspace(0,1,Ylims(i+1)+1);
[XI YI] = meshgrid(xx,yy);
D{end-i} = interp2(X,Y,D{end - i+1},XI,YI);

% Расчет f в точках X,Y
f = -B*exp(-A*((XI-xx(end)/2).^2+(YI-yy(end)/2).^2));
[M N] = size(D{end-i});

% Прогоняем невязку через метод Якоби (1 итерация)
for m = 2:M-1
for n = 2:N-1
D{end-i}(m,n) = 1/4*(D{end-i}(m-1,n)+ D{end-i}(m+1,n) + D{end-i}(m,n - 1) + D{end-i}(m,n+1)) - delta^2 /4 * f(m,n);
end
end

end
% Прибавляем получившуюся невязку к начальному решению
Fi_end = Fi+D{1};
% Рассчитываем максимальную относительную разницу между полученным решением и решением с предыдущей итерации
delta = ((Fi_end-Fi)./Fi_end);
maxdelta = max(delta(:))
Fi = Fi_end;

end


Вернуться к началу
 Профиль Отправить личное сообщение  
Ответить с цитатой  
СообщениеДобавлено: 22 дек 2013, 20:49 
Не в сети

Зарегистрирован: 13 сен 2011, 12:23
Сообщения: 46
Откуда: Тамбов
Прошу прощения за несамостоятельность своих решений, Даниил Анатольевич, но возник еще вопрос: "Из чего лучше исходить при выборе условия выхода из итерационного цикла?"
а) Из оперативной памяти компьютера;
б) Из того, насколько порядков нужно снизить погрешность, как было в лекциях (тут мы опять же упираемся в память компьютера);
в) Из того, что решения, найденные разными методами должны слабо отличаться друг от друга (трудно сделать);
г) Задать deltaMax = 1e-6, как мы обычно делаем и обосновать это тем, такой точности при относительной погрешности хватит наверняка;
д) Нет правильного ответа.


Вернуться к началу
 Профиль Отправить личное сообщение  
Ответить с цитатой  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 2 ] 

Часовой пояс: UTC + 3 часа


Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и гости: 1


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Перейти:  
cron
Создано на основе phpBB® Forum Software © phpBB Group
Русская поддержка phpBB