# Solving a linear system Ku=f using conjugate gradient method I implmented CG and PCG . Need help with visualisation of results

13 visualizaciones (últimos 30 días)
Lokesh el 21 de Nov. de 2023
Comentada: Lokesh el 21 de Nov. de 2023
In my code I implemented conjugate and preconjugate gradient method and visualised the result. I want to visualise the norm of the difference of the solution in every iteration with the final solution to see the convergence. Any help is appreciated. Below is my code:
clear all;
% very large matrix
dim = 2000;
vec = -1*ones(dim-1,1);
K = 4*eye(dim) + diag(vec,1)+diag(vec,-1);
K(1,1) = 100; % change Matrix to get an illconditioned matrix
% compute Condition matrix
cond(K)
%define righthandside
f = ones(dim,1);
%set maximal iteration and Toleranz
maxiter = 1000;
tol = 10e-6;
%set initial vector
u_0 = zeros(dim,1);
%run your implemented cg method
[u_cg, iter] = cg(K,f,maxiter,tol,u_0);%%% change code here !!!!!!!!!!!!!
u=u_0;
r=f-u;
d=r;
for k = 0:maxiter
alpha = (r.'*r) / (d.'*K*d); % Step size
u = u + alpha*d; % Update solution
r_new = r - alpha*K*d; % Update residual
% Correct search direction
beta = (r_new.'*r_new) / (r.'*r);
d = r_new + beta*d;
r = r_new; % Update residual
% Check for convergence
if norm(r) < Tol
break;
end
end
u_cg = Solution_using_Cg;
k = number_of_iterations;
%run your implemented pcg
[u_pcg, iter_pcg] = p_cg(K,f,maxiter,tol,u_0); %%% change code here !!!!!!!!!!!!!
u = u_0; % Initial guess
r = f - K*u; % Initial residual
C = diag(1./diag(K)); % Preconditioner
h = C*r; % Preconditioned residual
d = h; % Initial search direction
for k = 0:maxiter
alpha = (r.'*h) / (d.'*K*d); % Step size
u = u + alpha*d; % Update solution
r_new = r - alpha*K*d; % Update residual
h_new = C*r_new; % Preconditioned new residual
beta = (r_new.'*h_new) / (r.'*h);
d = h_new + beta*d; % Update search direction
r = r_new; % Update residual
h = h_new; % Update preconditioned residual
% Check for convergence
if norm(r) < Tol
break;
end
end
u_pcg;
iter_pcg;
%compare the results
norm(u_cg-u_pcg)
%use matlab pcg function and compare results
upcg = pcg(K,f);
norm(u_pcg-upcg)
%use matlab pcg function and compare results
u_lin = linsolve(K,f);
norm(u_pcg-u_lin)
##### 0 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

recent works el 21 de Nov. de 2023
visualize the norm of the difference of the solution in every iteration with the final solution to see the convergence:
clear all;
% very large matrix
dim = 2000;
vec = -1*ones(dim-1,1);
K = 4*eye(dim) + diag(vec,1)+diag(vec,-1);
K(1,1) = 100; % change Matrix to get an illconditioned matrix
% compute Condition matrix
cond(K)
%define righthandside
f = ones(dim,1);
%set maximal iteration and Toleranz
maxiter = 1000;
tol = 10e-6;
%set initial vector
u_0 = zeros(dim,1);
% run your implemented cg method
[u_cg, iter] = cg(K,f,maxiter,tol,u_0);
% compute norm of the difference of the solution in every iteration with the final solution
norm_diff_cg = zeros(iter, 1);
for i = 1:iter
norm_diff_cg(i) = norm(u_cg(1:i) - u_cg);
end
% run your implemented pcg
[u_pcg, iter_pcg] = p_cg(K,f,maxiter,tol,u_0);
% compute norm of the difference of the solution in every iteration with the final solution
norm_diff_pcg = zeros(iter_pcg, 1);
for i = 1:iter_pcg
norm_diff_pcg(i) = norm(u_pcg(1:i) - u_pcg);
end
% visualize the norm of the difference
figure;
semilogy(1:iter, norm_diff_cg);
xlabel('Iteration');
ylabel('Norm of the Difference');
title('Norm of the Difference for CG');
figure;
semilogy(1:iter_pcg, norm_diff_pcg);
xlabel('Iteration');
ylabel('Norm of the Difference');
title('Norm of the Difference for PCG');
This code will generate two plots, one for the norm of the difference of the solution in every iteration with the final solution for the CG method, and the other for the PCG method. You can compare the two plots to see how the convergence rate differs between the two methods.
##### 1 comentarioMostrar -1 comentarios más antiguosOcultar -1 comentarios más antiguos
Lokesh el 21 de Nov. de 2023
It is showing 2 empty plots

Iniciar sesión para comentar.

### Categorías

Más información sobre Logical en Help Center y File Exchange.

R2023b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by