% ode_2nd.m
% -u'' + a*u' + b*u = f, (0,T)
% u(0) = alp, u(T) = beta
% Scheme :
% A(i)*V(i-1) + B(i)*V(i) + C(i)*V(i+1) = F(i), i=2,3, ... ,N-1
% A(i) = -1/h^2-a(i)/(2*h), B(i) = 2/h^2+b(i), Ci = -1/h^2+a(i)/(2*h)
% Define coefficient functions and exact solution
% function ode_2nd(N)
a = @(t) cos(t);
b = @(t) exp(-t).*sin(t);
u = @(t) cos(t).*sin(t);
du = @(t) 1-2*sin(t).^2;
ddu = @(t) -4*cos(t).*sin(t);
f = @(t) -ddu(t)+a(t).*du(t)+b(t).*u(t);
% parameters
N = 100;
T = pi;
h = T/(N-1);
alp = 0; bet = 0; % boundary values
% Nodes
t = linspace(0,T,N); t = t';
% Define Coefficient matrix M = tridiag(A B C)
A = -1/h^2 - a(t)./(2*h); % N-vectors
B = 2/h^2 + b(t);
C = -1/h^2 + a(t)./(2*h);
M = diag(A(3:N-1),-1) + diag(B(2:N-1)) + diag(C(2:N-2), 1);
% Define F
F = f(t);
F(2) = F(2) - A(2)*alp;
F(N-1) = F(N-1) - C(N-1)*bet;
% Find approximate solution V
V = zeros(N,1); % initialize
V(1) = alp; V(N) = bet;
V(2:N-1) = M \ F(2:N-1);
% Exact solution U
U = u(t);
% Make a plot and find ||U-V||
subplot(121); plot(t,U, t, V,'.-'); title('Solution vs Approximation')
subplot(122); plot(t,U-V); title('The error plot')
Error = h*norm(U-V); % L^2-error
fprintf(' L2-Error = %12.4E \n', Error)
% The result :
L2-Error = 3.9696E-005
|
  |
| LAST UPDATE: 2009.11.05 - 16:03 |
|
메일보내기
이름검색
수치해석 lecture
matlab lecture fdm
