FiniteElement.m
classdef FiniteElement %FINITEELEMENT Summary of this class goes here
properties N%number of nodes
H% step in x direction
domain % 0 to 1
beta% hyperparameter
%initial conditions
alpha %analytic solution
exact f%f(x, t)
%Define boundary conditions
g0 % alpha0
g1 %alphaN
g0dash %alpha'0
g1dash %alpha'N
%stiffness matrices
A D end methods function obj = FiniteElement(N, beta) % Constructor initialize all variables
%Parameters Initialization
obj.N = N; obj.beta = beta; obj.H = 1 / N; % step in x direction obj.domain = linspace(0, 1, N + 1); obj.f = @(x, t) (1 - pi^2 * obj.beta ) * exp(t) * sin(pi * x); %Set initial conditions here
obj.alpha = sin(pi * obj.domain(2:end - 1)); %Set exact solution here
obj.exact = @(t) exp(t) * sin(pi*obj.domain); %Define boundary conditions obj.g0 = @(t) 0; obj.g1 = @(t) 0; obj.g0dash = @(t) 0; obj.g1dash = @(t) 0; obj = obj.matrices_req(); end function output = phi(obj, i, x) % basis functions see live script for better representation
if i == 0 if x >= 0 && x <= obj.H output = -(x - 1 * obj.H) / obj.H; else output = 0; end elseif i == obj.N if x >= (obj.N - 1) * obj.H && x <= obj.N * obj.H output = (x - (obj.N - 1) * obj.H) / obj.H; else output = 0; end else if (x >= (i - 1) * obj.H) && (x <= i * obj.H) output = (x - (i - 1) * obj.H) / obj.H; elseif x >= i * obj.H && x <= (i + 1) * obj.H output = -(x - (i + 1) * obj.H) / obj.H; else output = 0; end end end function output = phidash(obj, i, x) % derivative
%derivative of phi
if i == 0 if x >= 0 && x <= obj.H output = -1 / obj.H; else output = 0; end elseif i == obj.N if x >= (obj.N - 1) * obj.H && x <= obj.N * obj.H output = 1 / obj.H; else output = 0; end else if x >= (i - 1) * obj.H && x <= i * obj.H output = 1 / obj.H; elseif x >= i * obj.H && x <= (i + 1) * obj.H output = -1 / obj.H; else output = 0; end end end function output = const(obj, i, t) % constant in the matrix equation
output = obj.g0dash(t) * integral(@(x) obj.phi(i, x) * obj.phi(0, x), 0, 1) + ... obj.g1dash(t) * integral(@(x) obj.phi(i, x) * obj.phi(obj.N, x), 0, 1) ... -obj.beta * obj.g0(t) * integral(@(x) obj.phidash(i, x) * obj.phidash(0, x), 0, 1) + ... -obj.beta * obj.g1(t) * integral(@(x) obj.phidash(i, x) * obj.phidash(obj.N, x), 0, 1)... +obj.beta*obj.phi(i, 1)*obj.phidash(obj.N, 1)*obj.g1(t)-obj.beta*obj.phi(i, 0)*obj.phidash(0, 0)*obj.g0(t) ; end function obj = matrices_req(obj) % stiffness matrices
obj.A = zeros(obj.N - 1); obj.D = zeros(obj.N - 1); for i = 1:obj.N - 1 for j = 1:obj.N - 1 if ismember(abs(i - j), [0, 1]) obj.A(i, j) = integral(@(x) obj.phi(i, x) * obj.phi(j, x), 0, 1); obj.D( i, j) = integral(@(x) obj.phidash(i, x) * obj.phidash(j, x), 0, 1); end end end end function derivative_ = dydx(obj, t, y) % derivative used in ode45
C = zeros(obj.N - 1, 1); F = zeros(obj.N - 1, 1); for i = 1:obj.N - 1 C(i) = obj.const(i, t); F(i) = integral(@(x) obj.phi(i, x) * obj.f(x, t), 0, 1); end derivative_ = obj.A \ (-(-obj.beta * obj.D * y + C) + F); end function output = exactsoln(obj, times) % give exactsolution the same structure as output from ode45
output = zeros(length(times), length(obj.domain)); for i = 1:length(times) output(i,:) = obj.exact(times(i)) ; end end function [timestamps, ApproxTemp, ExactTemp] = temperature(obj, time_at) % main function to be called to get temperatures
[timestamps, ApproxTemp] = ode45(@(t, y) obj.dydx(t, y), time_at, obj.alpha); ApproxTemp = horzcat( arrayfun(@(x) obj.g0(x),timestamps),... ApproxTemp,... arrayfun(@(x) obj.g1(x), timestamps) ); ExactTemp = obj.exactsoln(timestamps); end end endrun.m
obj = FiniteElement(15, -1);[timestamp, A, E] = obj.temperature(linspace(.01, .1, 10)); x = obj.domain;y = timestamp;[xx, yy] = meshgrid(x, y);subplot(2, 1, 1);z = A;heatmap(x, y, z);colormap(flipud(hot));title("Approximate Solution");xlabel("The rod as x-axis");ylabel("Time");subplot(2, 1, 2);z = E;heatmap(x, y, z);colormap(flipud(hot));title("Exact Solution");xlabel("The rod as x-axis");ylabel("Time");Error:Error in FiniteElement/phi (line 76) if (x >= (i - 1) * obj.H) && (x <= i * obj.H)Error in FiniteElement>@(x)obj.phi(i,x)*obj.phi(j,x) (line 140) obj.A(i, j) = integral(@(x) obj.phi(i, x) * obj.phi(j, x), 0, 1);Error in integralCalc/iterateScalarValued (line 314) fx = FUN(t);Error in integralCalc/vadapt (line 132) [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);Error in integralCalc (line 75) [q,errbnd] = vadapt(@AtoBInvTransform,interval);Error in integral (line 88)Q = integralCalc(fun,a,b,opstruct);Error in FiniteElement/matrices_req (line 140) obj.A(i, j) = integral(@(x) obj.phi(i, x) * obj.phi(j, x), 0, 1);Error in FiniteElement (line 53) obj = obj.matrices_req(); >> ylabel("Time");
Best Answer