heat2.m

algrithm for optimization matlab source code for heat conducting plate 3 methods are included

function [t,x,U,Temp] = heat2 (T,a,m,n,beta,f,g1,g2)
%----------------------------------------------------------------------
% (d/dt)u(t,x) = beta*(d/dx)*(d/dx)u(t,x)
%
% The heat equation is to be solved over the region:
%
% 0 <= t <= T,
% 0 < x < a
%
% subject to the following initial and boundary conditions:
%
% u(0,x) = f(x)
% u(t,0) = g1(t)
% u(t,a) = g2(t)
%
%
% Inputs: T = final solution time (T > 0)
% a = upper boundary of x (a > 0)
% m = number steps in t (m >= 1)
% n = number of steps in x (n >= 1)
% beta = heat equation coefficient (beta > 0)
% f = string containing name of user-supplied
% initial condition function.
% g1,g2 = string containing name of user-supplied
% boundary condtion function
%
% The functions f and g are of the form:
%
% function y = f(x)
% function y = g(t)
%
%
% U = (m+1) by n matrix containing the estimated solution
% with U(k,j) = u(t(k),x(j)).
%
%----------------------------------------------------------------------
% Initialize
U = zeros (m+1,n);
Temp=zeros(m+1,n+2);
dt = T/m;
dx = a/(n+1);
t = [0 : dt : T]';
x = [0 : dx : (n+1)*dx]';
lambda = beta*dt/(dx*dx);
% Initial Condition
% u(0,x) = f(x)
for j = 1 : n
U(1,j) = feval(f,j*dx);
end
Temp(1,:)=[feval(g1,dt) U(1,:) feval(g2,dt)];
% Compute coefficient matrix of u(k+1)
A=zeros(n,n);
c0 = (1 + 2*lambda);
for i = 2 : n-1
A(i,i)=c0;
A(i,i-1)=-lambda;
A(i,i+1)=-lambda;
end
A(1,1)=c0;
A(n,n)=c0;
A(1,2)=-lambda;
A(n,n-1)=-lambda;
% Solve the PDE
for k = 1 : m
b = U(k,:);
s1=feval(g1,(k+1)*dt);
s2=feval(g2,(k+1)*dt);
b(1) = b(1) + lambda * s1;
b(n) = b(n) + lambda * s2;
U(k+1,:)=(A\b')';
Temp(k+1,:)=[s1 U(k+1,:) s2];
end
plot(x,Temp);
%----------------------------------------------------------------------
function [t,x,U,Temp] = heat3 (T,a,m,n,beta,f,g1,g2)
% Explicit
%----------------------------------------------------------------------
%----------------------------------------------------------------------
% Initialize
U = zeros (m+1,n);
Temp=zeros(m+1,n+2);
dt = T/m;
dx = a/(n+1);
t = [0 : dt : T]';
x = [0 : dx : (n+1)*dx]';
lambda = beta*dt/(dx*dx);
% Initial Condition
% u(0,x) = f(x)
for j = 1 : n
U(1,j) = feval(f,j*dx);
end
Temp(1,:)=[feval(g1,dt) U(1,:) feval(g2,dt)];
% Solve the PDE
c0=1-2*lambda;
for k = 1 : m
for j= 2 : n-1
U(k+1,j) = lambda * (U(k,j-1) + U(k,j+1)) + c0 * U(k,j);
end
U(k+1,1) = lambda * (feval(g1,k*dt) + U(k,j+1)) + c0 * U(k,j);
U(k+1,n) = lambda* (U(k,j-1)+feval(g2,(k+1)*dt))+ c0 * U(k,j);
Temp(k+1,:)=[feval(g1,(k+1)*dt) U(k+1,:) feval(g2,(k+1)*dt)];
end
plot(x,Temp);
%----------------------------------------------------------------------
function [t,x,U,Temp] = heat (T,a,m,n,beta,f,g1,g2)
%----------------------------------------------------------------------
% Crank-Nicolson
% Initialize
U = zeros (m+1,n);
Temp=zeros(m+1,n+2);
dt = T/m;
dx = a/(n+1);
t = [0 : dt : T]';
x = [0 : dx : (n+1)*dx]';
lambda = beta*dt/(dx*dx);
% Initial Condition
% u(0,x) = f(x)
for j = 1 : n
U(1,j) = feval(f,j*dx);
end
Temp(1,:)=[feval(g1,dt) U(1,:) feval(g2,dt)];
% Compute coefficient matrix of u(k+1)
A=zeros(n,n);
c0 = 2*(1 + lambda);
c1=2*(1-lambda);
for i = 2 : n-1
A(i,i)=c0;
A(i,i-1)=-lambda;
A(i,i+1)=-lambda;
end
A(1,1)=c0;
A(n,n)=c0;
A(1,2)=-lambda;
A(n,n-1)=-lambda;
% Solve the PDE
for k = 1 : m
for j = 2 : n-1
b(j) = lambda*(U(k,j-1) + U(k,j+1))+ c1*U(k,j);
end
s1=feval(g1,(k+1)*dt);
s2=feval(g2,(k+1)*dt);
b(1) = lambda*(feval(g1,k*dt) + U(k,2) + s1 )+ c1 * U(k,1);
b(n) = lambda*(U(k,n-1) + feval(g2,k*dt) + s2) + c1 * U(k,n);
U(k+1,:)=(A\b')';
Temp(k+1,:)=[s1 U(k+1,:) s2];
end
plot(x,Temp);
%----------------------------------------------------------------------