1-D Heat equation
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
I have no idea what went wrong. Two days since looking for the errors, please help
if
close all
x_min= 0;
x_max = 1;
N = 10;
L = x_max-x_min;
h = L/N;
x = linspace(x_min,x_max,N+1)'
uexact= @(x) 4*(x-x.^2);
f = @(x) 8*15; % Source term
u_exact = uexact(x);
ua = uexact(x_min);
ub = uexact(x_max);
u = zeros(N-1,1);
A = zeros(N-1,N-1);
b = f(x(2:N)); %
b(1) = b(1) + ua/h^2;
b(N-1) = b(N-1) + ub/h^2;
dA = diag( 2*ones(1,N-1) );
dAp1 = diag( -1*ones(1,N-2), 1 );
dAm1 = diag( -1*ones(1,N-2), -1 );
A = (dA + dAp1 + dAm1);
A = 15*A/h^2;
u = A\b; % solving the linear system
g = [ua; u; ub];
plot(x,u_exact,'b',x,g,'ro-');
legend('exact','numerical');
Thanks for the help.
0 commentaires
Réponse acceptée
Birdman
le 22 Mar 2018
x_min= 0;
x_max = 1;
N = 10;
L = x_max-x_min;
h = L/N;
x = linspace(x_min,x_max,N+1)'
uexact= @(x) 4*(x-x.^2);
f = @(x) 8*15*ones(1,N-1); % Source term
u_exact = uexact(x);
ua = uexact(x_min);
ub = uexact(x_max);
u = zeros(N-1,1);
A = zeros(N-1,N-1);
b = f(x(2:N)); %
b(1) = b(1) + ua/h^2;
for i=2:N-1
b(i) = b(i) + ub/h^2;
end
dA = diag( 2*ones(1,N-1) );
dAp1 = diag( -1*ones(1,N-2), 1 );
dAm1 = diag( -1*ones(1,N-2), -1 );
A = (dA + dAp1 + dAm1);
A = 15*A/h^2;
u = A\b.'; % solving the linear system
g = [ua; u; ub];
plot(x,u_exact,'b',x,g,'ro-');
legend('exact','numerical');
8 commentaires
Birdman
le 22 Mar 2018
This looks like a new question. Ask it to the forum as a new question so more people can contribute.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Digital Filtering dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!