Code not working.

1 次查看（过去 30 天）
adrooney
adrooney 2015-12-3
回答： Walter Roberson 2015-12-3
Hi, I got this code on the mit website, However I am not able to run the code please have a look
%---------------------------------------------------------------------
levels = 5; % size of problem
nu1 = 2; % number of presmoothing iterations
nu2 = 2; % number of postsmoothing iterations
gamma = 2; % number of coarse grid iterations (1=V-cycle, 2=W-cycle)
%---------------------------------------------------------------------
n = 2^(levels+2)-1; % number of grid points
h = 1/(n+1);
x = (h:h:(1-h))';
f = pi^2*(sin(pi*x)+4^2*sin(pi*4*x)+9^2*sin(pi*9*x));
A = spdiags(ones(n,1)*[-1 2 -1],-1:1,n,n);
b = f*h^2;
uc = A\b;
global t level, t = 0; level = levels;
clf, subplot(2,2,3:4), hold on
u = twogrid(A,b,nu1,nu2,gamma);
hold off, axis tight
subplot(2,2,1), plot(x,u,'b.-',x,uc,'r.--')
title('correct solution and multigrid approximation')
subplot(2,2,2), plot(x,uc-u,'r.-')
title('error')
%=====================================================================
function x = twogrid(A,b,nu1,nu2,gamma,x0)
%TWOGRID
% Recursive twogrid cycle for 1d Poisson problem.
% nu1 = number of presmoothing iterations (Gauss-Seidel)
% nu2 = number of postsmoothing iterations (Gauss-Seidel)
% gamma = number of coarse grid iterations (1=V-cycle, 2=W-cycle)
% x0 = starting vector (0 if not prescribed)
global t level
n = length(b);
if n<4
x = A\b; % solve exactly at coarsest level
else
G = speye(n)-tril(A)\A; cG = tril(A)\b; % create Gauss-Seidel
I = spdiags(ones(n-2,1)*[1 2 1],-2:0,n,n-2); % create interpolation
I = I(:,1:2:end)/2; R = I'/2; % and restriction matrices
if nargin<6, x = b*0; else, x = x0; end % starting vector
for i = 1:nu1, x = G*x+cG; end % presmoothing
r = b-A*x; % compute residual
rh = R*r; % restrict residual to coarse grid
t = t+1; level = level-1; plot([t-1 t],[level+1 level],'bo-')
eh = rh*0; % starting vector
for i = 1:gamma
eh = twogrid(R*A*I,rh,nu1,nu2,gamma,eh); % coarse grid iteration
end
e = I*eh; % interpolate error
t = t+1; level = level+1; plot([t-1 t],[level-1 level],'bo-')
x = x+e; % update solution
for i = 1:nu2, x = G*x+cG; end % postsmoothing
end

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2015-12-3
It appears to work fine for me.
Note: in order to use a "function", the code for the function must be written in to a .m file and that .m file must start with the word "function" (or sometimes "classdef"). It is not allowed to enter a function at the command prompt, and it is not allowed to have a "function" in a .m file that starts with executable code such as an assignment.
You can store from "function x = twogrid" onward in a file named twogrid.m . Or you can store everything in one .m file, provided that you start that .m file with
function TheNameOfTheFileGoesHere
For example I used
function test258671
and stored it all in test258671.m

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Transportation Engineering 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by