% New way to solve ode
% 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)
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
t = linspace(0,T,N); t = t'; % Nodes
A = -1/h^2 - a(t)./(2*h); B = 2/h^2 + b(t); C = -1/h^2 + a(t)./(2*h);
M = diag(A(2:N),-1) + diag(B(1:N)) + diag(C(1:N-1), 1);
M(1,:) = 0; M(1,1) = 1;
M(N,:) = 0; M(N,N) = 1;
F = f(t);
F(1) = alp; F(N) = bet;
% Find approximate solution V
V = M \ F;
% 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)
|
|
메일보내기
이름검색
ode45
수치해석 lecture
