math-bvp

% program FFT DEMO

% a script to walk you through the main features of the main features of
% matlab’s fft routines
%
close all; clear;
%
% Define an interval [a, b] on which the function is periodic, i.e.
% f(x + 2L) = f(x), where 2L = b-a
a = 0;
b = 2*pi;
twoL = b-a; L = twoL/2;
%
% Set the total number of grid points, M should be 2^p and initialize the
% Fourier Modes, k, and wave number
M = 16; N = M/2;
k = -N:N-1;
wk = pi*k/L;
%
% Construct x_j, make sure you don’t include the periodic point!
% x_1 = a
dx = (b-a)/M;
xj = a + 0:dx:b-dx;
%
% define a periodic function
f = inline(‘2 + sin(pi*x/L) + 0.5*cos(3*pi*x/L)’,’x’,’L’);
% f = inline(‘sin(pi*x/L)’,’x’,’L’);
fj = f(xj,L);
figure(1)
plot (a:.01:b,f(a:.01:b,L))
hold on
plot (xj,fj,’ro’)
title (‘Function f(x)’)
xlabel (‘x’)
legend (‘f(x)’,’f_j’,’Location’,’SouthWest’)
print -djpeg90 function.jpg
%
% Expected Fourier Coefficients
c0 = 2;
c1 = -i/2;
cm1 = i/2;
c3 = 1/4;
cm3 = 1/4;
%
% Calculate the Fourier coefficients
ck = exp(-i*wk*a).*fftshift(fft(fj))/M;
figure(2)
plot (k,real(ck),’ro’)
hold on
plot (k,imag(ck),’bo’)
plot (0,c0,’r*’)
hold on
plot (3,c3,’r*’)
plot (-3,cm3,’r*’)
plot (1,imag(c1),’b*’)
plot (-1,imag(cm1),’b*’)
legend (‘real(ck) – fft’,’imag(ck) – fft’,…
‘Location’,’NorthWest’)
title(‘Fourier coefficients’)
xlabel(‘k’)
print -djpeg90 fc_of_function.jpg
%
% Return to the original function using an ifft
fk = M*ifft(fftshift(exp(i*wk*a).*ck))
figure(3)
plot (xj,fk,’bo’)
hold on
plot (xj,fj,’r*’)
xlabel(‘x’)
ylabel(‘f(x)’)
title(‘Using ifft to return to original function values’)
print -djpeg90 ifft_of_fc.jpg

% Solve the two point boundary value problem using 2nd order accurate
% finite differences.
% u”(x) = f(x), x=[a,b]
% This was used to produce the plots shown in class on 27/09/12.
%
close all; clear;
%
% Define boundaries and rhs
f = inline(‘(2+4*x.^2).*exp(x.^2)’,’x’)
u_exact = inline(‘exp(x.^2)’,’x’)
a = 0; b = 1;
%
% Define an array of decreasing values of h
igrid = 5:10;
%
% for each value of h, set up grid on [0, 2Pi] and calculate D and f’
for j = igrid
m = 2^j;
h = (b-a)/m; xj = a:h:b;
hgrid(j-4) = h;
%
% evaluate f on grid;
fj = f(xj); % f_j = f(x_j)
fj(1) = u_exact(a); fj(end) = u_exact(b); % Dirichlet BCs
%
% Construct diagonal and offdiagonal entries for A
A = zeros(m+1,m+1);
diagonal = -2*ones(m+1,1)/h^2;
upper = (1/h^2)*ones(m,1); lower = (1/h^2)*ones(m,1);
A = diag(diagonal) + diag(upper,1) + diag(lower,-1); % tridiagonal matrix
A(1,1) = 1; A(1,2) = 0; % fix for left bc
A(m+1,m+1) = 1; A(m+1,m) = 0; % fix for right bc
%
% get norm of inverse
Ainv = inv(A);
normA(j-4) = norm(Ainv);
%
% solve
u = Afj’;
%
% plot numerical solution against exact
figure(1)
plot (xj, u_exact(xj),’r’)
hold on
plot (xj, u, ‘b’)
xlabel(‘x’)
ylabel(‘u’)
title (‘Solution to uprimeprime (x) = f(x)’)
legend (‘Exact’,’Computed’,’Location’,’NorthWest’)
if j == 5
print -djpeg90 ‘twopoint_bvp_dir.jpg’
end
pause; close all;
%
% calculate error
error(j-4) = norm(u’-u_exact(xj));
%
end
figure(2)
loglog(hgrid,error)
hold on
xlabel(‘h’)
ylabel(‘error’)
title(‘Error in Computed Solution’)
print -djpeg90 ‘err_2pntbvp_dir.jpg’
figure(3)
loglog(hgrid,normA)
xlabel(‘h’)
ylabel(‘||inv(A)||’)
print -djpeg90 ‘stability_2pntbvp_dir.jpg’
%
% Find solutions to the nonlinear BVP
% u” + u (u’-1) a < x < b
% u(a) = alpha u(b) = beta
%
close all; clear;
%
% Boundary conditions (Dirichlet)
alpha = 1; beta = 5;
%
% Define physical grid and calculate an initial guess
m = input(‘ Enter m: ‘);
a = 0; b = 1;
h = (b-a)/m; xj = a + h*(0:m);
u = alpha + (beta-alpha)*xj/(b-a);
plot (xj,u)
xlabel(‘ x’); ylabel(‘ u’);
title(‘ iterative solutions, u^k’)
hold on
%
% Define Newton Iteration Parameters
tol = 1.d-6; del = 1.d-4;
itmax = 100;
err = 1; iter = 0;
%
% Iteration loop
while err>tol
iter = iter + 1;
%
% Construct residual vector
F = zeros(m+1,1);
F(1) = 0;
F(2:m) = (u(3:m+1) – 2*u(2:m) + u(1:m-1))/h^2 …
+ u(2:m).*(u(3:m+1)-u(1:m-1))/(2*h) – u(2:m);
F(m+1) = 0;
%
% Assemble Jacobian matrix
jac = zeros(m+1,m+1);
diagonal = [1;-2*ones(m-1,1)/h^2 + (u(3:m+1)’-u(1:m-1)’)/(2*h) …
– ones(m-1,1);1];
upper = [0; ones(m-1,1)/h^2 + u(2:m)’/(2*h)];
lower = [ones(m-1,1)/h^2 – u(2:m)’/(2*h); 0];
jac = diag(diagonal) + diag(upper,1) + diag(lower,-1);
%
% If we are on the first iteration, check the Jacobian
if (iter==1)
for icol = 2:m
temp = u;
tempF = F;
delta_i = [zeros(icol-1,1);del;zeros(m+1-icol,1)];
u = u + delta_i’;
F = zeros(m+1,1);
F(1) = 0;
F(2:m) = (u(3:m+1) – 2*u(2:m) + u(1:m-1))/h^2 …
+ u(2:m).*(u(3:m+1)-u(1:m-1))/(2*h) – u(2:m);
F(m+1) = 0;
dF = (F-tempF)/del;
err_jac = norm(jac(:,icol) – dF);
disp([‘ column ‘ num2str(icol) …
‘ error in jacobian = ‘ num2str(err_jac)]);
u = temp;
F = tempF;
end
end
jac = sparse(jac);
delta = -jacF;
u = u + delta’;
err = norm(delta);
disp([‘it ‘ int2str(iter) ‘ error = ‘ num2str(err)])
plot (xj,u)
end;
if (iter>=itmax)
disp(‘ ITMAX EXCEEDED!’)
else
disp([‘CONVERGED IN ‘ int2str(iter) ‘ iterations’])
end;
plot (xj,u,’r’)

WE ARE A LEADING RESEARCH FIRM OFFERING AFFORDABLE RESEARCH SERVICES TO STUDENTS ON ALL SUBJECTS. ORDER NOW AND GET YOUR PAPER IN AS LITTLE AS 3 HOURS!!

You can place an order similar to this with us. You are assured of an authentic custom paper delivered within the given deadline besides our 24/7 customer support all through.