From: Peter Schaefer Date: Mon, 2 Apr 2012 18:18:01 +0000 (+0200) Subject: [src] erweitere A_plots auf mehere Dateien X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=7decb64a3311889062e12a772d5aff5b4bd0a958;p=bacc.git [src] erweitere A_plots auf mehere Dateien [src] Conditionierung Implementiert --- diff --git a/doc/doc.pdf b/doc/doc.pdf index cf94a25..fd4522c 100644 Binary files a/doc/doc.pdf and b/doc/doc.pdf differ diff --git a/doc/doc.tex b/doc/doc.tex index 1029755..e5497be 100644 --- 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 ae0219c..12188c3 100644 --- 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 0e208ce..6f3fa5c 100644 --- 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;