epsilon_small.m
8daa7c7f64acmaster

epsilon_small.m
View Options


			E = [1e-1, 1e-2, 1e-3]
			H = [1./16, 1./64, 1./512]
			alpha = 1
			beta = 3
			fconst = 1

			errors = zeros(3,3)

			for ie = 1:3
			for ih = 1:3
			epsilon = E(ie)
			h = H(ih)

			N = ceil(1.0/h);
			h = 1.0 / (N+1);

			% define matrix A and vector f
			A = zeros(N,N);
			f = ones(N,1) .* fconst;

			% compose matrix A
			for n = 1:N
			if n > 1
			A(n,n-1) = - (1.0/(2*h) + epsilon/h^2);
			end

			A(n,n) = 2 * epsilon / h^2;

			if n < N
			A(n,n+1) = 1.0/(2*h) - epsilon/h^2;
			end
			end

			% set boundary conditions
			f(1) = f(1) + (1.0/(2h) + epsilon/h^2) alpha;
			f(N) = f(N) - (1.0/(2h) - epsilon/h^2) beta;

			% solve system of linear equations
			u = A \ f;

			% compute analytical solution
			x = zeros(N,1);
			for i = 1:N
			x(i) = i * h;
			end

			ex = exact(alpha, beta, epsilon, x);

			% compare two results
			errors(ie,ih) = max(abs(u - ex))

			% output plots
			figure
			plot(x, u, '-*', x, ex, '-o')
			title(sprintf('Epsilon = %f, h = %f', epsilon, h))
			xlabel('x')
			ylabel('u')
			legend('Simulated','Exact')
			end
			end

			function ex = exact(alpha, beta, epsilon, x)
			ex = alpha + x + (beta - alpha - 1) * (exp( -(1 - x)/epsilon) - exp(-1.0/epsilon)) / (1 - exp(- 1.0/epsilon));
			end