% Bandmatrismetoden/Finita differensmetoden för ickelinjärt problem (GNM
% Exempel 6.24)

% Temperaturfördelning i glödtråd (symmetrisk kring x = 0.5)
% u'' = sigma*u^4 - I^2*(1 + 0.08*u), u(0) = 10, u'(0.5) = 0, sigma = 1e-7
% Bestäm u(x) för I = 75, 100, 125, 150 A

sigma = 1e-7;

figure(1)
clf
hold on
xlabel('X', 'FontSize', 14)
ylabel('Temperatur', 'FontSize', 14)


disp('  I   u(0.5)')

for I = [75 100 125 150]
	
	if (I == 75)
		col = 'r';
	elseif (I == 100)
		col = 'g';
	elseif (I == 125)
		col = 'b';
	else
		col = 'c';
	end
	
	
	% välj n punkter x_(i) = i*h med tillhörande värden u_(i)  
	% u_(0) = u(0), u_(n) = u_(0.5)
	% OBS! Tar med den ändpunkt som har ett derivatarandvillkor!
	n = 100;
	h = (0.5 - 0)/n;
	u0 = 10;
	% u'(0.5) = 0 => u_(n+1) = u_(n-1)
	
	% u'_(n+1) = 0  =>  u_(n) = u_(n + 2) 
	
	% Centraldifferens u_(i)'  = (u_(i+1) -  u_(i-1))/(2h)
	%                  u_(i)'' = (u_(i+1) - 2u_(i) + u_(i-1))/h^2
	% i = 1         : u0 - 2u_(1) + u_(2) = 
	%                   h^2(sigma*u_(1)^4 - I^2*(1 + 0.08*u_(1)))
	%                 (0.08*h^2*I^2 - 2)*u_(1) - sigma*h^2*u_(1)^4 + u_(2)
	%                   + (u0 + h^2*I^2) = 0
	% i = 2,..,(n-1): u_(i-1) - 2u_(i) + u_(i+1) = 
	%                   h^2(sigma*u_(i)^4 - I^2*(1 + 0.08*u_(i)))
	%                 u_(i-1) + (0.08*h^2*I^2 - 2)*u_(i) - sigma*h^2*u_(i)^4
	%                   + u_(i+1) + h^2*I^2 = 0
	% i = n         : u_(n-1) - 2u_(n) + u_(n-1) = 
	%                   h^2(sigma*u_(n)^4 - I^2*(1 + 0.08*u_(n)))
	%                 2u_(n-1) + (0.08*h^2*I^2 - 2)*u_(i) - sigma*h^2*u_(i)^4
	%                   + h^2*I^2 = 0
	% System: G(u) = Au + F(u) = 0, där A är en bandmatris (tridiagonal)
	A = spdiags([ones(n,1) (0.08*h^2*I^2 - 2)*ones(n,1) ones(n,1)], -1:1, n, n);
	A(n, n-1) = 2;

	% Startvärde för Newton-Raphson: Börja med att bestämma maxtemperaturen 
	% u(0.5) approximativt. Anta att temperaturen inte ändras i mitten av 
	% tråden, dvs u''(0.5) = u'(0.5) = 0
  % sigma*u05^4 - I^2*(1 + 0.08*u05) = 0
	u05 = fzero(@(u) sigma*u^4 - I^2*(1 + 0.08*u), 2000);
	% Anta sen att temperaturen växer från u0 till u05 för x = 0 till 0.1,
	% sen konstant
	ustart = [u0 + (u05 - u0)*(1:(n/5))/(n/5), u05*ones(1, 4*n/5)]';
	
	% Newton-Raphson
	Hnorm = n;
	u = ustart;
	while Hnorm > norm(u)*1e-8
		F = h^2*I^2 - sigma*h^2*u.^4;
		F(1,1) = F(1,1) + u0;
		G = A*u + F;
		% Jakobianen J = dG/du = A + dF/du
		% dF/du diagonal eftersom inga u' i den olinjära delen
    dFdu = spdiags(-4*sigma*h^2*u.^3, 0, n, n);
		J = A + dFdu;
		% Steget i gradientens riktning, flerdimensionella motsvarigheten till
		% steget f/f'
		H = J\G;
		Hnorm = norm(H);
		u = u - H;
	end
		
	fprintf(' %3.0f %8.3f\n', I, u(end));

	x = (0:h:1)';
	plot(x, [u0; u; flipud(u(1:(end - 1))); u0], col, 'LineWidth', 2)
	pause
	
end
	
legend('I = 75 A', 'I = 100 A', 'I = 125 A', 'I = 150 A', ...
	'Location', 'South') 
