I m getting this strange error and really I can't make any sense out of it. Please help. I was using simpsons() rule with 3*obj.N nodes but I saw that integration results were bad(where actual integration was 15 it was coming 14.867..) so I decided to use the inbuilt integral() method.
simpsons.m(used for integration for)
function I = simpsons(f,a,b,n)
h=(b-a)/n;
s = f(a) + f(b);
for i=1:2:n-1
s = s+ 4*f(a+i*h);
end
for i=2:2:n-2
s = s + 2*f(a+i*h);
end
I = h/3 * s;
end
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
end
run.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");
Aucun commentaire:
Enregistrer un commentaire