%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Syfte:    Löser, med cG(1)-metoden, differentialekvationen%
%%
%%          -(a(x)*u'(x))' + c(x)*u(x) = f(x),    x_min < x < x_max,%
%%
%%           med Robinrandvillkor%
%%
%%   a(x_min)*u'(x_min) = gamma(x_min)*(u(x_min) - g_D(x_min)) + g_N(x_min)
%%   -a(x_max)*u'(x_max) = gamma(x_max)*(u(x_max) - g_D(x_max)) + g_N(x_max)
%%
%%           och "plottar" lösningen.%
%%
%% Variabler:
%% p:  mesh           (Vektor med dimension 1 x N)
%% A:  styvhetsmatris (Matris med dimension N x N)
%% Mc:  massmatris    (Matris med dimension N x N)
%% R:  randmatris     (Matris med dimension N x N)
%% b:  lastvektor     (Vektor med dimension N x 1)
%% rv:  randvektor    (Vektor med dimension N x 1)
%% U:  lösning        (Vektor med dimension N x 1)
%%
%% Funktioner: a.m  definierar funktionen       a(x), x_min < x < x_max
%%             c.m  definierar funktionen       c(x), x_min < x < x_max
%%             f.m  definierar funktionen       f(x), x_min < x < x_max
%%
%%             gamma.m  definierar funktionen   gamma(x),   x = x_min, x_max
%%             g_D.m  definierar funktionen     g_D(x),   x = x_min, x_max
%%             g_N.m  definierar funktionen     g_N(x),   x = x_min, x_max

%%
%% MyFirstPoissonAssembler.m  beräknar A, Mc, R, b, och rv
%%
%%

%% Exempel:  Du vill lösa ekvationen
%%             -u''(x) + x*u(x) = (1 + x)*sin(x),
%%             med randvillkoren
%%             u(0) = 0,
%%             u'(pi/2) = 0,
%%             på en likformig indelning av [0, pi/2] i 9 delintervall.

%%             Börja med att definiera funktionerna a(x) = 1, c(x) = x,
%%             och f(x) = (1 + x)*sin(x), i filerna a.m, c.m, och f.m.

%%             Därefter definierar du funktionerna gamma, g_D, och g_N,
%%             i filerna gamma.m, g_D.m, och g_N.m.

%%             Tips: Dessa funktioner, som definierar randvillkoren, är endast
%%             definierade i randpunkterna x = x_min och x = x_max. 
%%             Inargumentet får alltså bara anta dessa två värden, och
%%             funktionerna skall returnera motsvarande funktionsvärde.
%%             Funktionen gamma skulle därför kunna tänkas konstrueras
%%             på följande sätt:

%%             function y = gamma(x)
%%             x_min = 0;
%%             x_max = pi/2;
%%             if x == x_min
%%              y = 10000; % Dirichletrandvillkor (approximativt)
%%             elseif x == x_max
%%              y = 0;     % Neumannrandvillkor
%%             else
%%              disp('Felaktigt argument i gamma.m')
%%             end

%% Det finns dock ett problem med konstruktionen ovan,
%% nämligen att det alltid är riskabelt att avgöra ifall det råder
%% likhet mellan två flyttal eftersom dessa inte representeras exakt i 
%% en dator. Ett annat (potentiellt) problem är att det skulle bli svårt att
%% generalisera rutinen till högre dimension, där randen är mer
%% komplicerad än bara två punkter. Ett bättre sätt (som också
%% enkelt generaliseras till högre dimension) är att ge varje del
%% av randen (i detta fall varje punkt) ett id-nummer:

%% function y = gamma(rand_nr)
%%   if rand_nr == 1       % vi ger x_min nummer 1
%%     y = 10000;       % (approximativt) Dirichletrandvillkor
%%   elseif rand_nr == 2   % vi ger x_max nummer 2
%%     y = 0;           % Neumannrandvillkor
%%   else
%%     disp('Felaktigt argument i gamma.m')
%%   end

%% Definiera nu på samma sätt g_D.m, och g_N.m, så att du får de önskade
%% randvillkoren. (Kommentar: Eftersom gamma(x_max) = 0, spelar det ingen
%% roll vilket värde du ger g_D(x_max).)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all        % Tar bl.a. bort alla variabler från Matlab:s workspace.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             %
% Definiera meshen p.         %
%                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


p = 0:5/9:5; %Notera att 9 delintervall motsvarar N = 10 nodpunkter
                         %(inkl. randnoderna).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                    %
% Assemblera matriser och vektorer.  %
%                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


[A, Mc, R, b, rv] = MyFirstPoissonAssembler(p);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             %
% Lös ekvationssystemet.      %
%                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


VL = A + Mc + R;
hl = b + rv;

U = VL\hl;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                    %
% Plotta lösningen.                  %
%                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  


plot(p, U)
xlabel('x')
title('U(x)')




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                           % 
% I detta fall är den exakta lösningen (verifiera detta!)   %
%                                                           %
%     u(x) = u = -e^x + x + e^5 - 5                         %
%                                                           %
% så du kan beräkna felet, u - U, enligt                    %
%                                                           %
% >>   u = -exp(p') + p' + exp(5) - 5;                      %
% >> fel = u - U;                                           %
%                                                           %
% (Vi transponerar 'p', eftersom 'U' är en kolonnvektor     %
%  och 'p' en radvektor.)                                   %
%                                                           %
% och plotta detta:                                         %
%                                                           %
% >> plot(p, fel)                                           %
%                                                           %
% Skriver du                                                %
%                                                           %
% >> max(abs(fel))                                          %
%                                                           %
% får du en uppskattning av felet i max-norm                %
% (d.v.s. max över alla noder).                             %
%                                                           %
% Det är mycket viktigt att testa sitt program genom att    %
% lösa några problem med kända (exakta) lösningar, och      %
% därefter jämföra den exakta med den beräknade lösningen.  %
%                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
