2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




 
 Старая тема
Сообщение12.12.2014, 13:04 
Аватара пользователя
Добрый день!

Тут было тема http://dxdy.ru/topic37436.html
Я слаб в методе прогонки чуток причесал, оно работает, но явно не правильно.
Помогите пожалуйста разобраться.
код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
clear, clc
L = 10; %Максимальное расстояние
T = 10; %Максимальное время
N=11; %Просто число...номер варианта=N
kappa=1;
%Устанавливаем шаг
h = L/(N-1); %по координате
tau = 0.01; %по времени
 
x = [0:h:L];
s = [0:tau:T];
%Граничные условия
Temp1=N*20;
Temp2=20*(N^2+N);
 
% Инициализируем сеточную функцию и обнуляем её
%Начальные условия
Temperature = zeros(length(s), length(x));
Temperature(1, :) = (x>=0).*(x<L/2)*Temp1 + (x>=L/2).*(x<=L)*Temp2;
Temperature(:, 1) = 220;
Temperature(:, N) = 341;
 
%Задаём коэффициенты для метода прогонки
A=kappa./h.^2;
C=1./tau-2*kappa./h.^2;
B=kappa./h.^2;
% Устанавливаем начальные прогоночные коэффициенты
F = (1/tau)*Temperature(2,2);
alpha(3) = -kappa./h.^2;
betta(3)=-Temp1*kappa./h.^2-(1./kappa)*Temperature(2,N);
 
%Запускаем обсчет по всем временам s
for i=1:length(s)
    for j=3:N-1 %прямая прогонка
        F = (1/tau)*Temperature(j,i);
        %Пригоночные коэффициенты
        alpha(j+1)=B/(C - alpha(j)*A);
        betta(j+1)=(A*betta(j) + F)/(C - alpha(j)*A);
    end
    Temperature(N,i+1) = Temp2; %Правая граница
    Temperature(N-1,i+1)=-((betta(N-1).*A+A*Temp2+(Temperature(j+1,i)./tau)))./((A*alpha(N-1)+2*A-1./tau));
    for j=N-1:-1:4, %обратная прогонка
        Temperature(j-2,i+1)=alpha(j-1)*Temperature(N-1,i+1)+betta(j-1);
    end;
end;
plot(s, Temperature(:,500),'-r') % Здесь надо T от времени нарисовать
 

 
 
 [ 1 сообщение ] 


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group