Matlab backup 4

clear;
clc;

% Parâmetros globais
eps = 10^-10;
Sro = 0+eps;      % Saturação residual do óleo
Srw = 0+eps;      % Saturação residual da água
phi = 1;        % Porosidade
V = 1;          % Velocidade
L = 3;          % Comprimento do domínio
T = 1;          % Tempo p

pasta = 'simulacoes/Linear';
mkdir(pasta);

% Função fw
fw = @(Sw) Sw;

% Derivadas das funções
dfw = @(Sw) 1;

% Parâmetros para a solução numérica
nx = 100;        % Número de pontos da malha espacial
nt = 100;        % Número de pontos da malha temporal
max_iter = 100;  % Máximo de iterações
tol = 1e-6;      % Tolerância

dx = L / (nx - 1);
x = linspace(dx/2, L-dx/2, nx-1);  % Malha espacial centralizada

dt = T/nt;
t = 0:dt:T;  % Malha temporal

% Inicialização da solução anterior com as condições iniciais
S_ant = Srw * ones(nx-1, 1);

% Adimensionalização
x_adim = x/L;
t_adim = (t*V) / (phi*L);

% Tempos que serão plotados
t_ref = (phi*L) / V;            % Tempo característico
T_estrela = T / t_ref;          % Tempo final adimensional

% Saturações
fileSaturacao = fopen([pasta '/saturacao.dat'], 'w');
fprintf(fileSaturacao,'%%');
fprintf(fileSaturacao,' Passo de tempo: %5d,  Tempo: %15.10f \n',0,0);
fprintf(fileSaturacao,'%25.20f %25.20f \n', [S_ant, S_ant]);
% Número de iterações
fileIteracoes = fopen([pasta '/Iteracoes.dat'], 'w');
fprintf(fileIteracoes,'%%');
fprintf(fileIteracoes,' Passo de tempo: %5d,  Tempo: %15.10f \n',0, 0);
fprintf(fileIteracoes,'%25.20f \n', 0);
% Número de saturação das iterações
fileSaturacaoIteracoes = fopen([pasta '/saturacaoIteracoes.dat'], 'w');
% Tempo
fileTempo = fopen([pasta '/tempo.dat'], 'w');
fprintf(fileTempo, '%25.20f \n', t(1));
% Espaço
fileEspaco = fopen([pasta '/espaco.dat'], 'w');
fprintf(fileEspaco, '%25.20f \n', x);
% Função de fluxo
fileFuncao = fopen([pasta '/funcao.dat'], 'w');
fprintf(fileFuncao, '%%');
fprintf(fileFuncao,' Passo de tempo: %5d,  Tempo: %15.10f \n',0,0);
fprintf(fileFuncao,'%25.20f %25.20f \n', [fw(S_ant), fw(S_ant)]);
% Derivada da função de fluxo
fileDerivada =  fopen([pasta '/derivada.dat'], 'w');
fprintf(fileDerivada, '%%');
fprintf(fileDerivada,' Passo de tempo: %5d,  Tempo: %15.10f \n',0,0);
fprintf(fileDerivada,'%25.20f %25.20f \n', [dfw(S_ant), dfw(S_ant)]);
% Linearização da função de fluxo
fileLinearizacao = fopen([pasta '/Linearizacao.dat'], 'w');

% Tempos que serão plotados
tempos = T_estrela * [0.5, 1];

% Método de Newton para a solução numérica
S_k1 = zeros(nx-1,1);
for j = 1:nt
    S_k = S_ant;

    fprintf(fileSaturacaoIteracoes,'%%');
    fprintf(fileSaturacaoIteracoes,' Passo de tempo: %5d,  Tempo: %15.10f \n',j,t(j+1));

    fprintf(fileLinearizacao,'%%');
    fprintf(fileLinearizacao,' Passo de tempo: %5d,  Tempo: %15.10f \n',j,t(j+1));

    for iter = 1:max_iter
        A = zeros(nx-1, nx-1);
        B = zeros(nx-1, 1);

        % Loop no espaço
        for i = 1:nx-1
            if i == 1
                % Preenche A(1,1)
                A(1, 1) = 1 + ((V * dt) / (phi * dx)) * (fw(S_k(i) + dfw(S_k(i)) * (S_k1(i) - S_k(i))) / S_k(i));
                % Preenche B(1)
                B(1) = S_ant(i) + ((V * dt) / (phi * dx)) * (fw(1 - Sro));
            else
                % Preenche o resto da matriz A
                fw_i_1 = fw(S_k(i - 1));
                fwi = fw(S_k(i));
                dw_i_1 = dfw(S_k(i - 1));
                dwi = dfw(S_k(i));

                % Coeficientes da matriz A
                A(i, i - 1) = -((V * dt) / (phi * dx)) * ((fw_i_1 + dw_i_1 * (S_k(i-1) - S_ant(i-1))) / S_k(i - 1));
                A(i, i) = 1 + ((V * dt) / (phi * dx)) * ((fwi + dwi * (S_k(i) - S_ant(i))) / S_k(i));

                % Preencher o resto do vetor B
                B(i) = S_ant(i);
            end
        end

        % Resolução do sistema linear
        S_k1 = A \ B;

        % Verificação do critério de convergência
        q = norm((S_k1 - S_k) / S_k1);
        if q < tol
            break;
        end

        h1 = fw(S_k(1)) + dfw(S_k(1)) * (S_k(1) - S_ant(1));
        h2 = fw(S_k(nx/2)) + dfw(S_k(nx/2)) * (S_k(nx/2) - S_ant(nx/2));
        h3 = fw(S_k(end)) + dfw(S_k(1)) * (S_k(end) - S_ant(end));

        fprintf(fileLinearizacao,'%25.20f %25.20f %25.20f \n', [h1, h2, h3]);

        fprintf(fileSaturacao,'%%');
        fprintf(fileSaturacao,' Passo de tempo: %5d,  Tempo: %15.10f \n',j, t(j+1));
        fprintf(fileSaturacao,'%25.20f %25.20f \n', [S_k, S_k1]);

         fprintf(fileFuncao,'%%');
         fprintf(fileFuncao,' Passo de tempo: %5d,  Tempo: %15.10f \n',j, t(j+1));
         fprintf(fileFuncao,'%25.20f %25.20f \n', [fw(S_k), fw(S_k1)]);

         fprintf(fileDerivada,'%%');
         fprintf(fileDerivada,' Passo de tempo: %5d,  Tempo: %15.10f \n',j, t(j+1));
         fprintf(fileDerivada,'%25.20f %25.20f \n', [dfw(S_k), dfw(S_k1)]);

        S_k = S_k1;

        fprintf(fileSaturacaoIteracoes,'%25.20f %25.20f %25.20f \n', [S_k1(1), S_k1(nx/2), S_k1(end)]);
    end

    fprintf(fileSaturacaoIteracoes,'%25.20f %25.20f %25.20f \n', [S_k1(1), S_k1(nx/2), S_k1(end)]);

    fprintf(fileLinearizacao,'%25.20f %25.20f %25.20f \n', [S_k1(1), S_k1(nx/2), S_k1(end)]);

    fprintf(fileIteracoes,'%%');
    fprintf(fileIteracoes,' Passo de tempo: %5d,  Tempo: %15.10f \n',j, t(j+1));
    fprintf(fileIteracoes,'%25.20f \n', iter);

    S_ant = S_k;

    fprintf(fileTempo, '%25.20f \n', t(j+1));
end