% This file used to be called Backward_Heat_SVD.m  Changed it to a forward solver.
% Galerkin discretization of the Fredholm integral equation for the forward heat equation
% g(x) = \int_0^pi k(x,s)f(s)ds where the kernel k(x,s) is defined by
% k(x,s) = \frac{2}{pi}\sum_{n=1}^{\infty}e^{-n^2}sin(nx)sin(ns)
% We will truncate the sum in k(x,s) as a sum from n=1 to N
%
% Jennifer Mueller and Samuli Siltanen, October 2012

clear all

% Compute data at time Tmax
Tmax = .4;

% Define the number of terms in the kernel 
N=10;
% Define the number of partition points of our interval [0,pi]
M=128
% Form the partition
x=zeros(M+1,1);
for k=1:M+1
  x(k) = (k-1)*pi/M; %  so x1=0, x2=pi/M, ... , x_(M+1)=pi
end
dx = x(2) - x(1)

A=zeros(M,M);
for ii = 1:M
  for jj=1:M
     for n=1:N
	 % From integrating sinnx * cos ns:
      A(ii,jj) = A(ii,jj) + exp(-(n^2)*Tmax)/n^2*(cos(n*x(ii))-cos(n*x(ii+1)))*(cos(n*x(jj))-cos(n*x(jj+1)));
    end
  end
end
A=2/(pi*dx)*A;

cond(A)
[U,D,V] = svd(A);

figure
semilogy(diag(D),'k.','MarkerSize',11)
axis([1 129 10^-25 10^0])
hold on
semilogy(diag(D),'k')
set(gca,'ytick',[10^-25,10^-20,10^-15,10^-5,10^0],'fontsize',15)
set(gca,'xtick',[1,64,128],'fontsize',15)
set(gca,'PlotBoxAspectRatio' , [2 1 1])
%title('Singular values of A')


%  Next set up the initial vector f1
alpha=2;
K=10;
for ii=2:M+1
  f1(ii-1,1) = K*sin(alpha*x(ii));  
end

uDat = A*f1;

%  Actual solution for u0 = 10 sin 2x is u(x,t) = 10*exp(-4t)*sin 2x
uAct = 10*exp(-4*Tmax)*sin(2*x(2:end));

% compute the relative sup norm error
sn = max(uAct)
supE=(norm(uDat - uAct,inf))/sn


% Show simulated data
figure
hold on
plot(x(2:end),uDat,'b')
plot(x(2:end),uAct,'--')
title('Data at time T (solid)')  
  