# Paulche/num_analysis_labs

1 parent 0ffa60a commit 6b61547c65a8195d7972f95e40f0cc390f99f689 committed Apr 7, 2013
Showing with 96 additions and 0 deletions.
1. +11 −0 lab9/action.m
2. +37 −0 lab9/compute_matrix.m
3. +20 −0 lab9/linput.m
4. +23 −0 lab9/solve_system.m
5. +5 −0 lab9/try_a.m
11 lab9/action.m
 @@ -0,0 +1,11 @@ +function [M F] = action() + [x p q k a b t n f] = linput; + + [a b c F M] = compute_matrix; + + Y = solve_system(a,b,c,F); + + plot_x = x(1):t:x(2); + + plot(plot_x,Y); +end
37 lab9/compute_matrix.m
 @@ -0,0 +1,37 @@ +function [a b c F M] = compute_matrix() + [x p q k a b t n f] = linput; + + a = []; + b = [k(1) - k(2)/t]; + c = [k(2)/t]; + + b(n+1) = k(3) - k(4)/t; + a(n+1) = - k(4) / t; + + for i = 1:(n-1) + i_e = i+1; + + q_v = q(x(1) + i * t); + p_v = p(x(1) + i * t); + + a(i_e) = t^-2 - p_v / (2*t); + b(i_e) = q_v - 2*t^-2; + c(i_e) = t^-2 + p_v / (2*t); + end + + for i = 1:(n+1) + F(i,1) = f(x(1) + (i-1)*t); + end + + M(1,1) = b(1); + M(1,2) = c(1); + + M(n+1,n) = a(n+1); + M(n+1,n+1) = b(n+1); + + for i = 2:n + M(i,i-1) = a(i-1); + M(i,i) = b(i-1); + M(i,i+1) = c(i-1); + end +end
20 lab9/linput.m
 @@ -0,0 +1,20 @@ +function [x p q k a b t n f] = linput() + % + % * y''(x) - (x + 3)^2 * y'(x) - 1/(x + 3)^2 * y(x) = 3, + % * x in [0,1] + % * y(0) - y'(0) = 4/3 + % * y(1) = 3/4 + % + n = 50; + x = [0 1]; + + + p = @(x) -(x + 3)^2; + q = @(x) -(x + 3)^-2; + f = @(x) 3; + + t = abs(x(1) - x(2)) / n; + a = 4/3; + b = 3/4; + k = [1 -1 1 0]; +end
23 lab9/solve_system.m
 @@ -0,0 +1,23 @@ +function Y = solve_system(a,b,c,F) + n = length(a) - 1; + + Ap = [-c(1)/b(1)]; + Bp = [F(1)/b(1)]; + + try_a(Ap(1)); + + for i = 2:n + Ap(i) = - c(i) / (a(i) * Ap(i-1) + b(i)); + + try_a(Ap(i)); + + Bp(i) = (F(i) - a(i) * Bp(i-1))/(a(i) * Ap(i-1) + b(i)); + end + + + Y(n+1) = (F(n+1) - a(n+1)*Bp(n))/(a(n+1)*Ap(n) + b(n+1)); + + for i = n:-1:1 + Y(i) = Ap(i) * Y(i+1) + Bp(i); + end +end
5 lab9/try_a.m
 @@ -0,0 +1,5 @@ +function try_a(value) + if abs(value) >= 1 + error('Ap > 1') + end +end