[src] erweitere A_plots auf mehere Dateien

[src] Conditionierung Implementiert

diff --git a/doc/doc.pdf b/doc/doc.pdf
index cf94a25ed416e8393a25886ff49677bed02db60e..fd4522c2e6a8418399c8991bb1f80936ebcd07f2 100644 (file)
Binary files a/doc/doc.pdf and b/doc/doc.pdf differ

diff --git a/doc/doc.tex b/doc/doc.tex
index 102975581c2495af0a5e5a3cb182a8c20bf2d07f..e5497be1700a22f9bd0654a6ee7f1f2e970b2fba 100644 (file)
--- a/doc/doc.tex
+++ b/doc/doc.tex
@@ -173,6 +173,18 @@ Damit ist $\phi_l$ die Galerkinapproximation an $\phi$
 \end{defi}
 
 
+\subsection{Vorkonditionieren}
+\begin{align}
+ V \phi_{h/2} &= b\\
+ D &= diag(V)\\
+ A &= D \cdot V \cdot D\\
+ c & = D\cdot b\\
+ A\cdot y &= c\\
+ \phi_{h/2} &= D \cdot y\\
+\end{align}
+
+
+
 \subsection{Netz}
 Sei $\T_h = \{T_1,T_2,\dots,T_N\}$ eine Triangulierung von $\Gamma$.
 \begin{itemize}

diff --git a/src/A_plots.m b/src/A_plots.m
index ae0219cddde9cfc1c659ee261c2bdfd6f0fd3fdf..12188c3736f7a54f4701474a6247473ac6e8c0ad 100644 (file)
--- a/src/A_plots.m
+++ b/src/A_plots.m
@@ -8,15 +8,16 @@ X = [];
 leg0 = {};
 leg1 = {};
 leg2 = {};
+leg3 = {};
 sym = {};
 
-rows = 8;
+rows = 9;
 
 
 for i = 1:length(files)
    
    load(files{i});
-   data
+%    data
       
    [m n] = size(data);
    step = round(n/rows);
@@ -34,7 +35,7 @@ for i = 1:length(files)
    end
    
    G_D = [G_D data(:,2:end)];
-   X = [X repmat(data(:,1),1,step*rows)];
+   X = [X repmat(data(:,1),1,step)];
    
    p1 = find(files{i}=='/',1,'last')+1;
    p2 = find(files{i}(p1:end)=='_',1);
@@ -56,6 +57,10 @@ for i = 1:length(files)
           ['min hmax/max hmax ' l0 l1{i}]...
           ['min hmin/hmax ' l0 l1{i}]...
           }';
+      
+      leg3 = {leg3{:}...
+          ['cond(A_{fine})' l0 l1{i}]...
+          }';
       sym = {sym{:} type2sym{data(1,[2+(i-1)*rows])}}';
    end
    
@@ -91,19 +96,19 @@ shift = shift+shift/10;
 % shift2 = shift2+shift2/10;
 % error*(eta(l)-shift2)/error(l)
 
-loglog(X(:,[2+(0:rows-5)]+rows*i),[G_D(:,2+rows*i) ...
+loglog(repmat(X(:,i+1),1,4),[G_D(:,2+rows*i) ...
     G_D(:,2+1+rows*i)...*(G_D(k,2)-shift)/G_D(k,3)...*G_D(1,2)/G_D(1,2+1+rows*i) ...
 ...    sqrt(abs(sol - G_D(:,2+2+rows*i)))...*G_D(1,2)/sqrt(abs(sol - G_D(1,2+2+rows*i)))...
     sqrt(abs(sol - G_D(:,2+2+rows*i)))...*(G_D(k,2)-shift)/G_D(k,3)...
     G_D(:,2+3+rows*i)...*G_D(1,2)/G_D(1,2+3+rows*i) ...
-    ],sym{i+1});
+    ],type2sym{i+1});
 hold on
 for i = 1:step-1
-loglog(X(:,[2+(0:rows-5)]+rows*i),[G_D(:,2+rows*i) ...
+loglog(repmat(X(:,i+1),1,4),[G_D(:,2+rows*i) ...
     G_D(:,2+1+rows*i)...*G_D(1,2)/G_D(1,2+1+rows*i) ...
     sqrt(abs(sol - G_D(:,2+2+rows*i)))...*G_D(1,2)/sqrt(abs(sol - G_D(1,2+2+rows*i)))...
     G_D(:,2+3+rows*i)...*G_D(1,2)/G_D(1,2+3+rows*i) ...
-    ],sym{i+1});
+    ],type2sym{i+1});
     
 end
 loglog(X(:,1),[7*X(:,1).^(-1/2),3*X(:,1).^(-1/4),2*X(:,1).^(-3/4)],'-.')
@@ -120,10 +125,10 @@ print('-r600','-depsc',[printt '_error.eps'])
  
 %% Plotte eNorm
 figure(5)
-loglog(X(:,2),G_D(:,2+2),sym{1});
+loglog(repmat(X(:,1),1,1),G_D(:,2+2),type2sym{1});
 hold on
 for i = 1:step-1
-    loglog(X(:,2+i*rows),G_D(:,2+2+i*rows),sym{i+1});
+    loglog(repmat(X(:,i+1),1,1),G_D(:,2+2+i*rows),type2sym{i+1});
 end
 loglog(X(:,1),repmat(sol,size(X,1),1),'r-.')
 hold off
@@ -139,18 +144,18 @@ print('-r600','-depsc',[printt '_norm.eps'])
 %% Plotte HMIN HMAX
 figure(6)
 i=0;
-loglog(X(:,(1:3)+rows*i),[...
+loglog(repmat(X(:,i+1),1,3),[...
     G_D(:,2+4+rows*i)...
     G_D(:,2+5+rows*i)...
     G_D(:,2+6+rows*i)...
-    ],sym{i+1});
+    ],type2sym{i+1});
 hold on
 for i = 1:step-1
-loglog(X(:,(1:3)+rows*i),[...
+loglog(repmat(X(:,i+1),1,3),[...
     G_D(:,2+4+rows*i)...
     G_D(:,2+5+rows*i)...
     G_D(:,2+6+rows*i)...
-    ],sym{i+1});
+    ],type2sym{i+1});
 end
 loglog(X(:,1),[7*X(:,1).^(-1/2),3*X(:,1).^(-1/4),2*X(:,1).^(-3/4)],'-.')
 hold off
@@ -163,6 +168,31 @@ legend({leg2{:} ...
      } ,'location','southwest','box', 'off');
 
 print('-r600','-depsc',[printt '_hminmax.eps'])
+
+
+%% Plotte cond
+figure(7)
+i=0;
+loglog(repmat(X(:,i+1),1,1),[...
+    G_D(:,2+7+rows*i)...
+    ],type2sym{i+1});
+hold on
+for i = 1:step-1
+loglog(repmat(X(:,i+1),1,1),[...
+    G_D(:,2+7+rows*i)...
+    ],type2sym{i+1});
 end
+loglog(X(:,1),[7*X(:,1).^(-1/2),3*X(:,1).^(-1/4),2*X(:,1).^(-3/4)],'-.')
+hold off
+
+title('hmin hmax')
+xlabel('Elemente');
+ylabel('Schaetzer');
+legend({leg3{:} ...
+      'N^{-1/2}'  'N^{-1/4}'  'N^{-3/4}'...
+     } ,'location','southwest','box', 'off');
+
+print('-r600','-depsc',[printt '_hminmax.eps'])
+
 
 end
\ No newline at end of file

diff --git a/src/A_step.m b/src/A_step.m
index 0e208ce15a45419fdda6b16f0726785b28188f74..6f3fa5c17b629be384115b5d03afe97d6dbfcb20 100644 (file)
--- a/src/A_step.m
+++ b/src/A_step.m
@@ -56,10 +56,10 @@ time = zeros(1,3);
   %Alle MatrixBrechenungsArten mit dem selben Netz berechnen
   for i = 1:length(type)
     %Matrix aufbauen -> MEX
-    A_fine = mex_build_V(coo_fine,ele_fine,zeta,type(i));
+    V_fine = mex_build_V(coo_fine,ele_fine,zeta,type(i));
     
     %Testet auf Fehlerhafte Einträge (NaN +/-Inf)
-    [r c] = find(isnan(A_fine)~=isinf(A_fine));
+    [r c] = find(isnan(V_fine)~=isinf(V_fine));
     if(~isempty(r))
         figure(9)
         plotShape(coo_fine,ele_fine(unique([r c]),:),'')
@@ -68,15 +68,32 @@ time = zeros(1,3);
     end
       
     %Lösung Berechnen
-    x_fine = A_fine\b_fine;
+    x_fine = V_fine\b_fine;
+    
+    %Vorkonditionierte Lösung!
+%     D = diag(diag(V_fine.^(-1/2)));
+%     
+%     A = D * V_fine * D;
+%     c = D*b_fine;
+%     y = A\c;
+%     x_fine = D*y;
     
     % \tilde \mu ( \Pi h -h + L_2 )
     tmu = hmin.* b .* sum((x_fine(f2s)'-repmat(sum(x_fine(f2s)',1)/4,4,1)).^2)' /4;
 
     
     %Fehlerschätzer 2 aufbauen
-    A = mex_build_V(G_C,G_E,zeta,type(i));
-    x = A\b;
+    V = mex_build_V(G_C,G_E,zeta,type(i));
+    
+    
+    x = V\b;
+    
+%     D = diag(diag(V.^(-1/2)));
+%     
+%     A = D * V * D;
+%     c = D*b;
+%     y = A\c;
+%     x = D*y;
     
     xo_fine(f2s) = repmat(x,1,4);
     xd_fine = xo_fine'-x_fine;
@@ -96,17 +113,17 @@ time = zeros(1,3);
     
     eta = abs(xe_fine-xe);
     
-   save_A_fine{i} = A_fine;
+   save_A_fine{i} = V_fine;
    save_x_fine{i} = x_fine;
    
-   save_A{i} = A;
+   save_A{i} = V;
    save_x{i} = x;
    
     
     
     
     data = [data type(i) sqrt(sum(tmu)) sqrt(eta) xe sqrt(sum(mu))...
-        min(hmin)/max(hmax) min(hmax)/max(hmax) min(hmin./hmax)];
+        min(hmin)/max(hmax) min(hmax)/max(hmax) min(hmin./hmax) cond(V_fine)];
   end
   time(2) = toc;