Introduction to Numerical Methods and MATLAB Programming for Engineers

% Plot the deflection y of a beam of length L under tension T
% and a uniform load w. Uses the exact analytic solution of the
% boundary value problem. For T=0, the solution breaks down.
% For T (e.g. T=.1) close to zero using the exact solution leads to loss
% of precision.
% Coefficients for a structural steel I-beam W12x22 in lb. in.
E = 29000000
I = 121
L = 120
% Load and tension
w = 10000
T = 10000
% calculated coefficients
EI = E*I
a = T/EI
b = w/(2*EI)
K = exp(sqrt(a)*L)
d5 = 2*b/(a^2)
d2 = - d5*(1/(1+K))
d1 = d2*K
d3 = b/a
d4 = - d3*L
% grid
xx = linspace(0,L,101);
% deflection
yy = d1*exp(-sqrt(a)*xx) + d2*exp(sqrt(a)*xx) + d3*xx.^2 + d4*xx + d5;
%plot the result
plot(xx,yy,'r')

% Finite elements program
% adapted from "Applied Numerical Anaylsis" by L. Faussett, Prentice Hall.
% Solves u_xx + u_yy = 0 
% with fixed boundary conditions 
% A list of vertices V and a list of triangles T must be in memory.
% x and y must be row vectors containing coordinates of the vertices.
% The following variable must also be in memory:
% n is the number of vertices
% m is the number of internal vertices
% p is the number of triangles
% z input values at the vertices, including boundary values.

% Plot the boundary conditions and initial values of nodes
figure
trimesh(T,x,y,z)

%Define basis function for each node
for j = 1:n                       % for each node (finite element)  
    for k = 1:p                   % for each triangle
        for i = 1:3               % for each vertex of the triangle
            aa(i,1) = 1;
            aa(i,2) = V(T(k,i),1);
            aa(i,3) = V(T(k,i),2);
            if T(k,i) == j        % if the vertex goes with the  element
                bb(i,1) = 1;      % set the value to one.
            else
                bb(i,1) = 0;      % otherwise set it to zero.
            end
        end
        abc = aa\bb;      % solve for coefficients of function j on T_k.
        A(j,k) = abc(1);
        B(j,k) = abc(2);
        C(j,k) = abc(3);
    end
end

% Calculate the areas of the triangles
H = ones(1,p);  
for i = 1:p
   H(i) = .5*abs(det([x(T(i,1)), x(T(i,2)), x(T(i,3)); y(T(i,1)), y(T(i,2)), y(T(i,3)); 1, 1, 1]));
end

% This section solves for the values at the interior vertices
for i = 1:m
    for j = 1:m
        s(1:p) = B(i,1:p).*B(j,1:p) + C(i,1:p).*C(j,1:p);
        AAA(i,j) = s*H';
    end
end
for i = 1:m
    for j = m+1:n
        k = j - m;
        s(1:p) = B(i,1:p).*B(j,1:p) + C(i,1:p).*C(j,1:p);
        G(i,k) = s*H';
    end
end
d(1:m) = -G*z(m+1:n)';
z(1:m) = AAA\d';


% Plot the solution
figure
trimesh(T,x,y,z)

%myheatdisk
R=5;      % radius of disk
n=100;    % number of subintervals
h=R/n;    % step size
hh=h^2;   % h^2 is used in each loop
steps=500;
d = .1;  % radiation coefficient
ub=10;    % temperature of surroundings
uR=10;    % temperature at the edge
ub4=ub^4; % ub^4 is used in each loop
r=-h/2:h:R+h/2;   % grid points
g=(r-5).^2;       % equation of source of heat
u=ones(1,n+2)*10; % sets initial temperatrue profile
u(n+2) = uR;      % sets temperature at edge
plot(r,u)         % plots initial temperature profile
hold on; 
for j=1:steps    
    for i=2:n+1  
     u(i)=-hh*d*(u(i)^4-ub4)/3 + hh*g(i)/3 + h*(u(i+1)-u(i-1))/(6*r(i)) + (u(i-1)+u(i)+u(i+1))/3;
    end
    u(1)=u(2); % makes u'(0) = 0
    plot(r,u)  % plots temperature profiles
end
hold off;

function I = mylowerleft(f,a,b,c,d,m,n)
    % Computes an approximate integral of a function of two variables f(x,y).
    % Uses the lower-left corner to evaluate.
    % Inputs:
    %    a,b -- define the interval in x, namely a<x<b.
    %    c,d -- define the interval in y, namely c<y<d.
    %    m -- the number of (evenly spaced) intervals to use in x.
    %    n -- the number of (evenly spaced) intervals to use in y.
    % Output: The approximate value of the integral.
    %
    % Note: This program intentionally has a bug. Find and fix it.
    h = (b-a)/m;                            % step size in x direction
    k = (d-c)/n;                            % step size in y direction
    [X,Y] = meshgrid(a:h:b,c:k:d);          % sets up a grid
    Z=f(X,Y);
    mesh(X,Y,Z)                             % plot, just for fun
    S = sum(sum(Z(1:m,1:n)));
    I = S*h*k;
end

function [T , Y] = mymodeuler(f,tspan,y0,n)
    % Solves dy/dt = f(t,y) with initial condition y(a) = y0
    % on the interval [a,b] using n steps of the modified Euler's method.
    % Inputs: f -- a function f(t,y) that returns a column vector of the same
    %              length as y
    %         tspan -- a vector [a,b] with the start and end times
    %         y0 -- a column vector of the initial values, y(a) = y0
    %         n  -- number of steps to use
    % Outputs: T -- a n+1 column vector containing the times
    %          Y -- a (n+1) by d matrix where d is the length of y
    %               Y(j,i) gives the ith component of y at time T(j) 

    a = tspan(1); b = tspan(2); % parse starting and ending points
    h = (b-a)/n;        % step size
    t = a; T = a; % t is the current time and T will record all times
    y = y0;       % y is the current variable values, as a column vector
    Y = y0';      % Y will record the values at all steps, each in a row
    for i = 1:n
        k1 = h*f(t,y);      % y increment based on current time
        k2 = h*f(t+h,y+k1); % y increment based on next time
        y = y + .5*(k1+k2); % update y using the average
        t = a + i*h;        % The next time.
        T = [T; t];         % Record t and y into T and Y.
        Y = [Y; y'];
    end
end

% mypoisson
% Poisson on a Rectangle
% Solves u_xx + u_yy = f(x,y)
% with boundary conditions:
% u(a,y) = u(x,c) = 1
% u(b,y) = 1 + 5/8 y
% u(x,d) = 1 + 1/2 x
% with f(x,y) = .5 sin(x)sin(y)

% Define the rectangle and nodes
a = 0;
b = 10;
c = 0;
d = 8;
m = 15;
n = 10;
h = (b-a)/m;
k = (d-c)/n;
x = a:h:b;
y = c:k:d;

% Assign number of iterations and calculate constants
maxit = 0;
c1 = h^2/(2*(h^2 + k^2));
c2 = k^2/(2*(h^2 + k^2));
c3 = k^2*c1;

% Assign initial values and boundary conditions
u = zeros(m+1,n+1);
u(1,:) = 1;
u(:,1) = 1;
u(:,n+1) = 1 + 5*sin(x(:)*pi/b/2);
u(m+1,:) = 1 + (5/8)*y(:);

% Assign values of f(x,y) in the interior
f = zeros(m+1,n+1);
for i = 2:m
f(i,:) = .5*sin(x(i)).*sin(y(:));
end

% Iterate values of u in the interior
for j = 1:maxit
u(2:m,2:n) = c1*(u(2:m,3:n+1)+u(2:m,1:n-1)) + c2*(u(3:m+1,2:n)+u(1:m-1,2:n)) - c3*f(2:m,2:n);
end

% Plot the result
mesh(x,y,u')

% Program to produce the triangulation,
% boundary condition and initial values
% at the interior nodes for a Washer.

clear all

% V contains vertices
V = [% interior points
    1.4712    0.2926
    1.2472    0.8334
    0.8334    1.2472
    0.2926    1.4712
   -0.2926    1.4712
   -0.8334    1.2472
   -1.2472    0.8334
   -1.4712    0.2926
   -1.4712   -0.2926
   -1.2472   -0.8334
   -0.8334   -1.2472
   -0.2926   -1.4712
    0.2926   -1.4712
    0.8334   -1.2472
    2.0000         0
    1.8478    0.7654
    1.4142    1.4142
    0.7654    1.8478
    0.0000    2.0000
   -0.7654    1.8478
   -1.4142    1.4142
   -1.8478    0.7654
   -2.0000    0.0000
   -1.8478   -0.7654
   -1.4142   -1.4142
   -0.7654   -1.8478
   -0.0000   -2.0000
    0.7654   -1.8478
    1.4142   -1.4142
    1.8478   -0.7654
    1.2472   -0.8334
    1.4712   -0.2926
    % outerboundary
    2.4520    0.4877
    2.0787    1.3889
    1.3889    2.0787
    0.4877    2.4520
   -0.4877    2.4520
   -1.3889    2.0787
   -2.0787    1.3889
   -2.4520    0.4877
   -2.4520   -0.4877
   -2.0787   -1.3889
   -1.3889   -2.0787
   -0.4877   -2.4520
    0.4877   -2.4520
    1.3889   -2.0787
    2.0787   -1.3889
    2.4520   -0.4877
    % inner boundary
    1.0000      0
    0.9239    0.3827
    0.7071    0.7071
    0.3827    0.9239
    0.0000    1.0000
   -0.3827    0.9239
   -0.7071    0.7071
   -0.9239    0.3827
   -1.0000    0.0000
   -0.9239   -0.3827
   -0.7071   -0.7071
   -0.3827   -0.9239
   -0.0000   -1.0000
    0.3827   -0.9239
    0.7071   -0.7071
    0.9239   -0.3827];
   
% T assigns the traingles
T = [43    42    25
    11    10    59
    11    10    25
    23    40    41
    22    40    39
    22    23    40
    21    38    39
    21    22     7
    21    22    39
    24    10     9
    24    42    41
    24    42    25
    24    10    25
    24    23    41
    24    23     9
    27    45    44
    28    45    46
    28    27    45
    30    47    48
    29    47    46
    29    28    14
    29    28    46
    29    30    47
    26    11    12
    26    43    44
    26    43    25
    26    11    25
    26    27    44
    26    27    12
    20    38    37
    20    21    38
    58    57     9
    58    10    59
    58    10     9
     8    57    56
     8    22    23
     8    23     9
     8    57     9
     8    56     7
     8    22     7
    63    62    14
    60    61    12
    60    11    59
    60    11    12
    13    61    62
    13    28    27
    13    27    12
    13    61    12
    13    62    14
    13    28    14
    19    36    37
    19    20    37
    19    20     5
    17    34    35
    55    56     7
    53    54     5
    18    36    35
    18    19    36
    18    17    35
    16    33    34
    16    17    34
     6    55    54
     6    20    21
     6    21     7
     6    55     7
     6    20     5
     6    54     5
     4    19     5
     4    53     5
     4    53    52
     4    18    19
     3    18    17
     3     4    52
     3     4    18
     2    16    17
     2     3    17
    31    29    30
    31    32    30
    31    32    64
    31    64    63
    31    29    14
    31    63    14
    15    30    48
    15    32    30
    15    33    48
    15    16    33
     1     2    16
     1    15    16
     1    15    32
    51     3    52
    51     2     3
    50     1     2
    50    51     2
    49     1    32
    49    32    64
    49    50     1];


% x, y are row vectors containing coordinates of vertices
x = V(:,1)';
y = V(:,2)';
% n is the number of vertices
n= size(V,1);
% m is the number of interior vertices
m = 32;
% p is the number of triangles
p= size(T,1);


% Plot the mesh
figure
trimesh(T,x,y)

% Set values at nodes
% The first 32 entries are interior vertices, the final 32 are
% boundary vertices.

z = x+2.5;    % Set node values to a linear funtion.
z(49:64) = 0; % Set inner boundary nodes to zero.
z(1:32) = 0;  % Set interior nodes to zero.

% Plot the mesh with node values specified.
figure
trimesh(T,x,y,z)

% Construct a Wedge-shaped region and its triangulation

clear all

% Set up internal nodes (manually):
V = [ .89*cos(pi/8) .89*sin(pi/8);
	 1.9*cos(pi/8) 1.9*sin(pi/8);
	 3*cos(pi/12) 3*sin(pi/12);
	 3*cos(pi/6) 3*sin(pi/6) ;
	 4*cos(pi/16) 4*sin(pi/16);
	 4*cos(pi/8) 4*sin(pi/8);
	 4*cos(3*pi/16) 4*sin(3*pi/16) 
	 4.6*cos(.03*pi) 4.6*sin(.028*pi);
	 4.6*cos(.22*pi) 4.6*sin(.222*pi) ];
m= size(V,1); % m = number internal nodes

% Set up boundary nodes (manually) and append:
a = cos(pi/4);
V = [ V;  0 0 ; 1 0 ; 2 0 ; 3 0 ; 4 0 ; 5 0 ;
       a a ; 2*a 2*a; 3*a 3*a; 4*a 4*a; 5*a 5*a]; 
for i = 1:4
  x = 5*cos(pi*i/20);
  y = 5*sin(pi*i/20);
  V = [ V; x y];
end
n = size(V,1);  % n = total number nodes

% extract x and y coordinates of all nodes
x = V(:,1)';
y = V(:,2)';
% make triangulation
T = delaunay(x,y); 
p = size(T,1); % p = number of triangles

% Plot the mesh
figure
trimesh(T,x,y)

% Set node values to a function
z = sin(x/2)+(y/3).^2;
z(1:9) = 0;

% Plot the mesh with node values specified.
figure
trimesh(T,x,y,z)