How to solve for 1D non homogenous ODE by Finite element method

조회 수: 4 (최근 30일)
MUTHULINGAM
MUTHULINGAM 2024년 9월 10일
편집: Torsten 2024년 9월 10일
I am unable to input the dirichlet condition into this non homogenous ode. Please point out my mistake.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% -u" + u = - 5x + 4 0 < x <= 1 %%%%
%%%% u(0) = 0, u(1) = 2 %%%%
%%%% Exact solution u = (exp(x)*(3*exp(1) + 4))/(exp(2) - 1) - 5*x - (exp(-x)*(3*exp(1) + 4*exp(2)))/(exp(2) - 1) + 4 %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
format long;
N = 20;
x0 = 0;
xN = 1;
h = 1/N;
for j = 1:N+1
X(j) = x0 + (j-1)*h;
end
f = @(x)(- 5*x +4);
A = zeros(N+1,N+1);
F = zeros(N+1,1);
%%%% Local stiffness matrix %%%%
a = [1/h -1/h ; -1/h 1/h] ;
for i=1:N
phi1 = @(x)(X(i+1)-x)/h; %%%% linear basis function %%%%
phi2 = @(x)(x-X(i))/h; %%%% linear basis function %%%%
f1 = @(x)f(x)*phi1(x); %%%% integrand for load vector %%%%
f2 = @(x)f(x)*phi2(x); %%%% integrand for load vector %%%%
v(1,i) = gauss(f1,X(i),X(i+1),1); %%%% element wise values of
v(2,i) = gauss(f2,X(i),X(i+1),1); %%%% load vector
end
Unrecognized function or variable 'gauss'.
%%%% Assembling %%%%
for i=1:N
A([i i+1],[i i+1]) = A([i i+1],[i i+1]) + a;
F([i i+1],1) = F([i i+1],1) + v([1 2],i);
end
%%%% Dirichlet Boundary condition %%%%
% F(N+1,1) = F(N+1,1)+2;
fullnodes = [1:N+1];
%%%%% Dirichlet boundary condition %%%%%
freenodes=setdiff(fullnodes,[1]);
Uh = zeros(N+1,1);
Uh(N+1,1) = 2;
%%%% Approximate solution %%%%
Uh(freenodes)=A(freenodes,freenodes)\F(freenodes,1);
%%%% Exact solution %%%%
U = zeros(N+1,1);
for i =1:N+1
U(i) = (exp(X(i))*(3*exp(1) + 4))/(exp(2) - 1) - 5*X(i) - (exp(-X(i))*(3*exp(1) + 4*exp(2)))/(exp(2) - 1) + 4;
end
error = max(abs(U-Uh));
[U Uh]
plot(X,U,X,Uh,'o')

답변 (1개)

Torsten
Torsten 2024년 9월 10일
편집: Torsten 2024년 9월 10일
I usually set
A(1,1) = 1, F(1) = 0, A(N+1,N+1) = 1, F(N+1) = 2
and only loop from 2 to N here:
%%%% Assembling %%%%
for i=1:N
A([i i+1],[i i+1]) = A([i i+1],[i i+1]) + a;
F([i i+1],1) = F([i i+1],1) + v([1 2],i);
end
I don't know if this is "Finite-Element-like".

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

제품

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by