Numerical soluton to heat equation MATLAB

Post by **answerhappygod** » Wed Mar 30, 2022 9:20 am

: Numerical Soluton To Heat Equation Matlab 1 (88.34 KiB) Viewed 29 times

: Numerical Soluton To Heat Equation Matlab 2 (48.61 KiB) Viewed 29 times

: Numerical Soluton To Heat Equation Matlab 3 (69.1 KiB) Viewed 29 times

This week's problem will test your comparison between the numerical and analytical solutions to the heat equation. The problem parameters for the person in the arctic tundra problem which were used in Week 4 should be kept the same. During the tutorial, you should have developed a numerical solution to the heat equation using the Forward in Time Central in Space (FTCS) method. As a reminder, the algorithm can be described as below: Tņ+1 = o7i+1+(1 – 20)T? +oT-1 where: 0 cdt dr2 Note that c is the heat equation constant, and dr is the radius step size. Write a MATLAB function that takes as input nr, dt and timeend, which represent the number of points in the radius vector, time step size in seconds, and the ending time, respectively. Your function should check whether the inputs meet the Von Neumann stability requirement and output two variables - the flag for stability stabilityFlag ('1' for stable or 'O' for unstable) and the maximum absolute percentage error Error between the numerical solution above and the analytical solution derived last week (using Fourier series) at the timepoint closest to timeend (e.g. if dt=5 and timeend = 21, you would compute the error at t=20 s). If timeend is exactly half-way between two timepoints, you should compute the error at the later timepoint.
Function nr 1 function (stabilityFlag, Error] = tutorial5(nr, dt, timeend) 2 % Description: This file provides the solution within which you can 3 % implement the analytical vs numerical solution to the heat equation 4 % 5 % Input parameters: 6 % - number of grid points in r 7 % dt - time step size in seconds 6 - timeend - ending time 9 % 10 % Output parameters: 11 % - stabilityFlag - the Boolean for stability 12 % - Error - maximum absolute percentage error of numerical vs analytical 13 %% YOU NEED TO FILL IN THE LINES WITH QUOTATION MARKS 00 14 = 15 % Spatial (x) parameters 16 L = 0.12; % Size of the domain 17 rstart - 0; % Start of computational domain (m) 18 rend = rstart + L; % End of computational domain (m) 19 dr (rend - rstart)/(nr - 1); % Spatial step size, delta r (m) 20 r = rstart: dr : rend; % Vector of grid points 21 22 % Phyiscal parameters 23 1 = -41; % Temperature at r=0 24 T2 = 37; % Temperature at r=L 25 A T2-T1; % Amplitude of the saw tooth wave 26 k = 0.527; % Thermal conductivity (W/m-K) 27 rho 1000; % Density (kg/m^3) 28 gamma = 3600; % Heat capacity (3/kg-K) 29 C = sqrt(k/rho/gamma); % Heat equation constant 의 = 30 31 % Temporal (t) parameters 32 timestart 0; % Starting time ( seconds) 33 time = timestart: dt : timeend; % vector for time (seconds) 34
35 % Initialise the solution arrays - Tn (Tan) and Tnp1 (T^(n+1)) and exact 36 % solution Texact 37 Tn zeros(1, nr); % Current value of temperature at time n 38 Tnp1 zeros (1, nr); % Next value of temperature at time n+1 39 TExact = zeros(1, nr); % Exact fourier series solution for temperature = = 40 = 41 % Specify the initial conditions 42 for j=2:nr-1 43 Inj) = ''; % initial conditions 44 end 45 46 % Enforce the boundary conditions - These will never change in the 47 % simulation 48 Tn(1) 49 Tn(nr) 50 Tnp1(1) 51 Tnp1(nr) 52 53 % Fourier parameters pre-compute B and lambda vectors 54 nfs = 1000; % Number of terms in the Fourier Series Expansion 55 B = zeros(1, nfs); 56 lambda = zeros(1, nfs); 57 for n = 1:nfs; B(n)=''; % Calculate B for each fourier term lambda (n)=""; % Calculate lambda for each fourier term 60 end = 58 59
61 62 %% Main solution loop 63 64 % Loop through time 65 for i-2:nt = % loop starts from 2 (not 1) to ensure that the analytical and numerical solutions are compared at the same time point in the mean error calc below 66 67 t=time(i); % current time for outputted solution 68 69 % Analytical solution (same as Week 3) 70 for j=1:nr; % for each grid point 71 TExact(j)=''; % add the steady state solution for n-1: nfs; % sum each individual term in the series 73 TExact(); 74 end end 76 77 % Numerical solution - note that we only need to solve for j=2 to nr-1 78 % since the first and last points are fixed by the boundary condition sigma-''; for j-2:nr-1 81 Tnp1(j)=''; 82 end 72 75 79 80 83 84 end 85 Code to call your function 1 nr = 51; 2 dt = 0.06; 3 timeend = 10; 4 [stabilityFlag, Error] = tutorial5(nr, dt, timeend); () =