From: Nils Forssén Date: Wed, 12 Oct 2022 14:05:17 +0000 (+0200) Subject: push X-Git-Url: https://gitweb.forssennils.se/?a=commitdiff_plain;h=ccb1229b61cba9340f29eb4690b2e1d00f2d6f45;p=TANA21.git push --- diff --git a/LorenzAttractor.m b/LorenzAttractor.m new file mode 100644 index 0000000..bc7ee53 --- /dev/null +++ b/LorenzAttractor.m @@ -0,0 +1,63 @@ +function LorenzAttractor(T) +%% +% Nonlinear ODE example from nonlinear dynamics to solve the +% Lorenz system using a Runge-Kutta 4 time integrator +% +% d(y1)/dt = 10*(y2-y1) +% d(y2)/dt = 28*y1 - y2 - y1*y3 +% d(y3)/dt = y1*y2 - 8/3*By3 +% +% INPUT: +% +% T - final time +% +% OUTPUT: +% +% NONE - this makes a movie +% +%% +% adhoc tactic to set the time step + num_steps = 2500; + dt = T/num_steps; +%% +% set the initial condition and allocate memory for the solution vector + y = zeros(3,num_steps); + ini = 5*ones(3,1); +%% +% save the initial condition and enter the Runge-Kutta loop to solve the ODE + y(:,1) = ini; + for i = 2:num_steps + k1 = RHS(y(:,i-1)); + k2 = RHS(y(:,i-1) + 0.5*dt.*k1); + k3 = RHS(y(:,i-1) + 0.5*dt.*k2); + k4 = RHS(y(:,i-1) + dt.*k3); + y(:,i) = y(:,i-1) + (dt/6).*(k1 + 2.*k2 + 2.*k3 + k4); + end +%% +% visualize the solution progessively to make a movie + clf; + hold on; + xlim([-20 20]); + ylim([-30 30]); + zlim([0 50]); + grid on; + for n = 1:(length(y)/7) + plot3(y(1,1:7*n),y(2,1:7*n),y(3,1:7*n),'b'); + drawnow + view(-37.5-n, 24) + pause(0.02) + end +end + +%%% +% auxiliary function to get the right hand side +%%% + function dy = RHS(y) +%% +% evaluate the right hand side of the nonlinear ODE for the Lorenz +% attractor. The input y is a vector with three components. + dy = zeros(3,1); + dy(1) = + dy(2) = + dy(3) = + end \ No newline at end of file diff --git a/Project/N_20.fig b/Project/N_20.fig new file mode 100644 index 0000000..5eb85ad Binary files /dev/null and b/Project/N_20.fig differ diff --git a/Project/N_20.svg b/Project/N_20.svg new file mode 100644 index 0000000..64b7199 --- /dev/null +++ b/Project/N_20.svg @@ -0,0 +1,180 @@ + + +11.11.21.31.41.51.61.71.81.920.50.520.540.560.580.60.620.640.660.680.7 diff --git a/Project/N_5.fig b/Project/N_5.fig new file mode 100644 index 0000000..d622a05 Binary files /dev/null and b/Project/N_5.fig differ diff --git a/Project/N_5.svg b/Project/N_5.svg new file mode 100644 index 0000000..f5e0e04 --- /dev/null +++ b/Project/N_5.svg @@ -0,0 +1,180 @@ + + +11.11.21.31.41.51.61.71.81.920.50.520.540.560.580.60.620.640.660.680.7 diff --git a/Project/lufact.m b/Project/lufact.m new file mode 100644 index 0000000..5ca3886 --- /dev/null +++ b/Project/lufact.m @@ -0,0 +1,49 @@ +function [L,U,P] = lufact(A) +%% +% Computes the LU factorization of the matrix A. The factorization is done +% in-place and then the L and U matrices are extracted at the end for +% output. The factorization is PA = LU +% +% INPUT: +% A - n by n square, non-singular matrix +% +% OUTPUT: +% L - lower triangular matrix with ones along the main diagonal +% U - upper triangular matrix +% P - permutation matrix for pivoting +%% +% get the size of the system + n = length(A); +%% +% initialize the pivoting matrix + P = eye(n); +%% +% Vectorized in-place LU factorization (with row pivoting) that keeps +% track of the total permutations by scrambleing the matrix P + for j = 1:n-1 + A + [v i] = max(abs(A(j:n, j))); + A([j i+j-1],:) = A([i+j-1 j],:); + P([j i+j-1],:) = P([i+j-1 j],:); + A +%% +% Find the index of the largest pivot element in the current column + + +%% +% Swap the rows within the in-place array as well as the permutation matrix P + + + +%% +% Perform the in-place elimination and save the new column of L + i = j+1:n; % indices for the "active" matrix portion + A(i,j) = A(i,j)/A(j,j); + A(i,i) = A(i,i) - A(i,j)*A(j,i); + %A, return + end +%% +% Extract L and U from the in-place form + U = triu(A); + L = eye(n) + tril(A,-1); +end diff --git a/Project/script.m b/Project/script.m new file mode 100644 index 0000000..1790f1f --- /dev/null +++ b/Project/script.m @@ -0,0 +1,59 @@ +%% Section 1 + +N = 5; + +% create three random vectors +main = N*ones(N,1) + rand(N,1); +upper = rand(N-1,1); +lower = rand(N-1,1); +% use the ’diag’ command to create tridiagonal matrix +A = diag(main,0) + diag(upper,1) - diag(lower,-1); + +known_x = ones(N,1); % make known solution vector of 1s +manuf_b = A*known_x; % manufacture the right-hand-side + +A + +%% Section 2 + +N = 5; +x = linspace(1,2, N); +h = x(2) - x(1); + +% Build tridiag matrix A +main = ones(N-2,1); +for i=2:N-1 + main(i-1,1) = (2/(x(i)^2)) - (2/(h^2)); +end +lower = ones(N-3,1); +for i=2:N-2 + lower(i-1,1) = (1/(h^2)) - (2/(x(i+1)*h)); +end +upper = ones(N-3,1); +for i=2:N-2 + upper(i-1,1) = (1/(h^2)) + (2/(x(i)*h)); +end + +A = diag(main,0) + diag(upper,1) + diag(lower,-1); + +% Build constant matrix b +b = zeros(N-2, 1); +for i=2:N-1 + b(i-1) = (2*log(x(i)))/(x(i)^2); +end +b(1) = b(1) + -((1/(h^2)) - (2/(x(1)*h)))*0.5; +b(end) = b(end) + -((1/(h^2)) + (2/(x(end)*h)))*log(2); + +answer = thomas(A, b); +answer = [0.5;answer;log(2)]; + +plot(x, answer); +hold on; +real_values = zeros(N, 1); +for i=1:N + real_values(i) = (4/x(i)) - (2/(x(i)^2)) + log(x(i)) - (3/2); +end +plot(x, real_values); + + +norm(answer - real_values, inf) \ No newline at end of file diff --git a/Project/thomas.asv b/Project/thomas.asv new file mode 100644 index 0000000..f4c60ae --- /dev/null +++ b/Project/thomas.asv @@ -0,0 +1,21 @@ +function [x] = thomas(A,b) + % from wikipedia + n = length(b); + c_prime = zeros(n-1, 1); + d_prime = zeros(n, 1); + c_prime(1) = A(1,2) / A(1,1); + d_prime(1) = b(1) / A(1,1); + for i=2:n-1 + c_prime(i) = A(i,1 + i) / (A(i,i) - (A(i, i-1)*c_prime(i-1))); + d_prime(i) = (b(i) - (A(i, i-1)*d_prime(i-1))) / (A(i, i) - (A(i, i-1)*c_prime(i-1))); + end + d_prime(n) = (b(n) - (A(n, n-1)*d_prime(n-1))) / (A(n, n) - (A(n, n-1)*c_prime(n-1))); + + % backward sub + x = zeros(n, 1); + x(n) = d_prime(n); + for j=n-1:-1:1 + x(j) = d_prime(j) - (c_prime(j)*x(j+1)); + end + +end \ No newline at end of file diff --git a/Project/thomas.m b/Project/thomas.m new file mode 100644 index 0000000..f4c60ae --- /dev/null +++ b/Project/thomas.m @@ -0,0 +1,21 @@ +function [x] = thomas(A,b) + % from wikipedia + n = length(b); + c_prime = zeros(n-1, 1); + d_prime = zeros(n, 1); + c_prime(1) = A(1,2) / A(1,1); + d_prime(1) = b(1) / A(1,1); + for i=2:n-1 + c_prime(i) = A(i,1 + i) / (A(i,i) - (A(i, i-1)*c_prime(i-1))); + d_prime(i) = (b(i) - (A(i, i-1)*d_prime(i-1))) / (A(i, i) - (A(i, i-1)*c_prime(i-1))); + end + d_prime(n) = (b(n) - (A(n, n-1)*d_prime(n-1))) / (A(n, n) - (A(n, n-1)*c_prime(n-1))); + + % backward sub + x = zeros(n, 1); + x(n) = d_prime(n); + for j=n-1:-1:1 + x(j) = d_prime(j) - (c_prime(j)*x(j+1)); + end + +end \ No newline at end of file diff --git a/RungeKutta4.m b/RungeKutta4.m new file mode 100644 index 0000000..c7cecff --- /dev/null +++ b/RungeKutta4.m @@ -0,0 +1,44 @@ +function [t,y] = RungeKutta4(f,t0,T,y0,h) +%% +% Classical four-stage Runge-Kutta time integration technique +% to solve the initial value problem +% +% y' = f(t,y), y(t0) = y0 +% +% INPUT: +% +% f - right-hand-side function that can depend on +% the time t and the function y +% t0 - intial time +% T - final time +% y0 - initial solution value +% h - fixed time step size +% +% OUTPUT: +% +% t - vector of time values +% y - solution vector of approximate values +% +%% +% discretize the time variable and set the problem size + t = t0:h:T; + y = zeros(length(t),1); +%% +% save the intial values + y(1) = y0; + t(1) = t0; +%% +% use Runge-Kutta 4 stage with the given timestep h to solve the problem and +% save the information for later analysis + for j = 2:length(t) + k_1 = f(t(j-1), y(j-1)); + k_2 = f(t(j-1) + (h/2), y(j-1) + ((h/2)*k_1)); + k_3 = f(t(j-1) + (h/2), y(j-1) + ((h/2)*k_2)); + k_4 = f(t(j-1) + h, y(j-1) + ((h*k_3))); + y(j) = y(j-1) + ((h/6)*(k_1 + (2*k_2) + (2*k_3) + k_4)); + + + + + end +end diff --git a/SimpsonsRule.m b/SimpsonsRule.m new file mode 100644 index 0000000..cecf9ad --- /dev/null +++ b/SimpsonsRule.m @@ -0,0 +1,72 @@ +function [S,x] = SimpsonsRule(f,a,b,n) +%% +% Implementation of the (possibly) composite Simpsons's quadrature +% rule to approximate the definite integral +% +% b +% /\ +% | +% | f(x) dx +% | +% \/ +% a +% +% INPUT: +% +% f - anonymous function integrand +% a - lower integration bound +% b - upper integration bound +% n - number of intervals (must be an even integer) +% +% OUTPUT: +% +% In - approximation to the definite integral on the given interval [a,b] +% x - set of n+1 quadrature nodes +% +%% +% Simpson's rule needs at least two intervals + if n == 1 + error('Simpsons rule requires n >= 2'); + end +%% +% check if the number of intervals is even + if mod(n,2) ~= 0 + error('The number of intervals n must be even'); + end +%% +% Find the interval spacing h + h = (b-a)/n; + +%% +% build the set of uniformly spaced quadrature nodes + x = zeros(n+1,1); + x(1) = a; + for i=2:n+1 + x(i) = x(i-1) + h; + +%% +% Initialize the integral value to zero + S = f(x(1)) + f(x(end)); + for i=2:n + if mod(i, 2) == 0 + S = S + 4*f(x(i)); + else + S = S + 2*f(x(i)); + end + end + + S = (h/3) * S; + +%% +% Approximate the integral with (composite) Simpson's rule + + + + + + + + + + +end diff --git a/adaptativeSimpson.m b/adaptativeSimpson.m new file mode 100644 index 0000000..5b96887 --- /dev/null +++ b/adaptativeSimpson.m @@ -0,0 +1,67 @@ +function [In,t] = adaptativeSimpson(f,a,b,tol) +%% +% Adaptive Simpson's rule to approximate a definite integral. This uses +% a posterioi error analysis to recursively apply Simpson's rule where it +% is needed +% +% INPUT: +% +% f - integrand function +% a - lower bound of interval +% b - upper bound of the interval +% tol - error tolerance for adaptivity +% +% OUTPUT: +% +% In - adaptive approximation of the definite integral to +% within the error tolerance +% t - vector of the adapted quadrature nodes +% +%% +%% +% use recursive bisection with error estimation to compute the integral approximation + m = 0.5*(b+a); % find the current midpoint + [In,t] = do_integral(a,f(a),b,f(b),m,f(m),tol); +%% +% this is a recursive function within Matlab so it calls itself + function [In,t] = do_integral(a,fa,b,fb,m,fm,tol) + %% + % need the two midpoints of the sub-intervals and the function + % evalutions for the recursion + xL = 0.5*(a+m); + fxL = f(xL); + xR = 0.5*(b+m); + fxR = f(xR); + %% + % save the five nodes at the current level of the recursion + t = [a;xL;m;xR;b]; + %% + % get the h value for the whole interval + h = 0.5*(b-a); + %% + % compute the Simpson's rule on the whole interval + S_coarse = h/3*(fa + 4*fm + fb); + %% + % compute the Simpson's rule on the two subintervals + S_Left = h/6*(fa + 4*fxL + fm); + S_Right = h/6*(fm + 4*fxR + fb); + S_fine = S_Left + S_Right; + %% + % error estimate + E = S_coarse - S_fine; + %% + % check against the user set tolerance + if abs(E) < 10*tol + % exit when tolerence is met + In = S_fine; + else + %% + % bisect again if the tolernce is note met + [IL,tL] = do_integral(a,fa,m,fm,xL,fxL,tol); + [IR,tR] = do_integral(m,fm,b,fb,xR,fxR,tol); + In = IL + IR; + % append the node edges together (this avoids duplicates) + t = [tL;tR(2:end)]; + end + end +end diff --git a/forwardEuler.m b/forwardEuler.m new file mode 100644 index 0000000..fa91267 --- /dev/null +++ b/forwardEuler.m @@ -0,0 +1,37 @@ +function [t,y] = forwardEuler(f,t0,T,y0,h) +%% +% Forward Euler time integration technique to solve the +% initial value problem +% +% y' = f(t,y), y(t0) = y0 +% +% INPUT: +% +% f - right-hand-side function that can depend on +% the time t and the function y +% t0 - intial time +% T - final time +% y0 - initial solution value +% h - fixed time step size +% +% OUTPUT: +% +% t - vector of time values +% y - solution vector of approximate values +% +%% +% discretize the time variable and set the problem size + t = t0:h:T; + y = zeros(length(t),1); +%% +% save the intial values + y(1) = y0; + t(1) = t0; +%% +% use forward Euler with the given timestep h to solve the problem and +% save the information for later analysis + for j = 2:length(t) +% this is an explicit method so we always use values from the current solution + y(j) = y(j-1) + h*f(t(j-1),y(j-1)); + end +end diff --git a/lab3_k.m b/lab3_k.m new file mode 100644 index 0000000..372f02f --- /dev/null +++ b/lab3_k.m @@ -0,0 +1,30 @@ +%lab3 +N = 1000; +x = linspace(0,1,N+1); +h = x(2) - x(1); +f_prime = zeros(1, N+1); + +f_1 = @(x) (1); +f_2 = @(x) (x); +f_3 = @(x) (x.^2); +f_4 = @(x) (exp(sin(4*x))); + +for i=1:N-1 + f_prime(i+1) = (f_4(x(i+2)) - f_4(x(i))) / (2*h); +end + +f_prime(1) = (f_4(x(2)) - (f_4(x(1)))) / (h); +f_prime(end) = (f_4(x(end)) - f_4(x(end-1))) / h; + +%f_prime(1) = ((-3*f_4(x(1))) + (4*f_4(x(2))) - (f_4(x(3)))) / (2*h); +%f_prime(end) = ((f_4(x(end-2))) - (4*f_4(x(end-1))) + (3*f_4(x(end)))) / (2*h); + +f_prime_exact = zeros(1, N+1); +for i=1:N+1 + f_prime_exact(i) = 4*exp(sin(4*x(i))) * cos(4*x(i)); +end + +f_prime +f_prime_exact + +max(abs(f_prime_exact - f_prime)) \ No newline at end of file diff --git a/plot_composite_quad.m b/plot_composite_quad.m new file mode 100644 index 0000000..d854aaa --- /dev/null +++ b/plot_composite_quad.m @@ -0,0 +1,24 @@ +function plot_composite_quad(f,t) +%% +% Script that will plot a function and the interval boundaries for use +% with a compositie quadrature rule +% +% INPUT: +% +% f - integrand function +% t - set of nodes used for the composite quadrature rule +% +% OUTPUT: +% +% NONE - Produces plot to screen +% +%% + clf + p1 = fplot(f,[t(1) t(end)],'-k','LineWidth',2); + hold on + p2 = stem(t,f(t),'fill','color',[0.9100 0.4100 0.1700],'MarkerSize',9,'LineWidth',1.5); + set(gcf,'Position',[500, 60, 1250, 1250]) + set(gca,'FontSize',16); + legend([p1 p2],{'$(x+1)^2\cos\left(\frac{2x+1}{x-3.3}\right)$','Quadrature nodes'}... + ,'interpreter','latex','fontsize',24,'Location','northwest') +end \ No newline at end of file