44 lines
1.2 KiB
Matlab
44 lines
1.2 KiB
Matlab
% Parameters
|
|
N = 501; % number of grid points
|
|
x = linspace(0,1,N); % grid
|
|
dx = x(2)-x(1); % grid spacing
|
|
p_values = [-10, -3, -1, 0, 1, 3, 10];
|
|
|
|
% Analytical solution function
|
|
u_analytical = @(x,p) (abs(p) < 1e-12).*x + (abs(p) >= 1e-12).*((exp(p.*x)-1)./(exp(p)-1));
|
|
|
|
figure;
|
|
hold on;
|
|
|
|
for k = 1:length(p_values)
|
|
p = p_values(k);
|
|
|
|
e = ones(N-2,1);
|
|
% Second derivative: -u''
|
|
D2 = spdiags([e -2*e e], -1:1, N-2, N-2)/dx^2;
|
|
% First derivative: p*u'
|
|
D1 = spdiags([-e*0.5 e*0.5], [-1 1], N-2, N-2)/dx;
|
|
L = -D2 + p*D1;
|
|
|
|
% Right-hand side
|
|
b = zeros(N-2,1);
|
|
|
|
% Incorporate boundary conditions
|
|
b(1) = b(1) + 0/dx^2; % u(0)=0
|
|
b(end) = b(end) + 1/dx^2; % u(1)=1
|
|
|
|
% Solve linear system
|
|
u_inner = L\b;
|
|
u_num = [0; u_inner; 1]; % add boundary values
|
|
|
|
% Plot numerical and analytical solution
|
|
plot(x, u_num, '--', 'LineWidth', 2, 'DisplayName', sprintf('num, p=%g',p));
|
|
plot(x, u_analytical(x,p), '-', 'LineWidth', 2, 'DisplayName', sprintf('ana, p=%g',p));
|
|
end
|
|
|
|
xlabel('x');
|
|
ylabel('u(x)');
|
|
title('Numerical vs Analytical solution of -u'''' + p u'' = 0');
|
|
legend show;
|
|
grid on;
|
|
hold off;
|