Chapter 28 The Solution with the 2nd-order RK (Heun without Corrector) Can Be Laid out As

CHAPTER 28 28.1 The solution with the 2nd-order RK (Heun without corrector) can be laid out as For the 4th-order RK, the solution is A plot of both solutions along with the analytical result is displayed below: 50 Analytical 25 RK-4 RK-2 0 0 20 40 28.2 The mass-balance equations can be written as dc1 = −0.14 c1 + 0.04c 3 dt dc 2 dt dc 3 dt dc 4 dt dc 5 dt = 0.2c1 − 0.2 c2 = 0.025c2 − 0.275c3 = 0.1125c3 − 0.175c4 + 0.025c5 = 0.03c1 + 0.03c2 − 0.06c5 Selected solution results (Euler’s method) are displayed below, along with a plot of the results. 20 10 0 0 20 c1 40 c2 c3 60 c4 80 c5 Finally, MATLAB can be used to determine the eigenvalues and eigenvectors: >> a=[.14 -.04 0 0 0;-.2 .2 0 0 0;0 -.025 .275 0 0;0 0 -.1125 .175 -.025;-.03 -.03 0 0 .06] a = 0.1400 -0.2000 0 0 -0.0300 -0.0400 0.2000 -0.0250 0 -0.0300 0 0 0.2750 -0.1125 0 0 0 0 0.1750 0 0 0 0 -0.0250 0.0600 >> [v,d]=eig(a) v = d = 0 0 0 1.0000 0 0 0 0.6644 -0.7474 0 0 0 0 0.2124 0.9772 -0.1836 -0.2954 -0.0370 0.1890 0.9176 0.0826 -0.2567 -0.6021 0.7510 0.0256 0.1750 0 0 0 0 0 0.2750 0 0 0 0 0 0.0600 0 0 0 0 0 0.0757 0 0 0 0 0 0.2643 28.3 Substituting the parameters into the differential equation gives dc = 20 − 0.1c − 0.1c 2 dt The mid-point method can be applied with the result: The results are approaching a value of 13.6351 Challenge question: The steady state form (i.e., dc/dt = 0) of the equation is 0 = 200 − c − c 2 , which can be solved for 13.65097141 and -14.65097141. Thus, there is a negative root. If we put in the initial y value as −14.650971 (or higher precision, the solution will stay at the negative root. However, if we pick a value that is slightly higher (a per machine precision), it will gravitate towards the positive root. For example if we use −14.65097 15 10 5 0 -5 0 -10 -15 -20 5 10 Conversely, if we use a slightly lower value, it will go unstable 0 0 5 10 -10 -20 -30 -40 28.4 The first steps of the solution are shown below along with a plot. Notice that the a value of the inflow concentration at the end of the interval (cin-end) is required to calculate the k2’s correctly. 50 0 0 50 100 28.5 The system is as depicted below: qevap = 0.5 kg/min qin = 10 kg/min sin = 8 g/kg qout = 10 kg/min M0 = 1000 kg s0 = 8 g/kg (a) The mass of water in the tank can be modeled with a simple mass balance dM = qin − qout − q evap = 10 − 10 − 0.5 = −0.5 dt With the initial condition that M = 1000 at t = 0, this equation can be integrated to yield, M = 1000 − 0.5t Thus, the time to empty the tank (M = 0) can be calculated as t = 1000/0.5 = 2000 minutes. (b) The concentration in the tank over this period can be computed in several ways. The simplest is to compute the mass of salt in the tank over time by solving the following differential equation: dm = qin sin − q out s dt where m = the mass of salt in the tank. The salt concentration in the tank, s, is the ratio of the mass of salt to the mass of water s= m m = M 1000 − 0.5t The first few steps of the solution of this ODE with Euler’s method is tabulated below. In addition, a graph of the entire solution is also displayed. 8.6 8.4 8.2 8 0 500 1000 1500 2000 Recognize that a singularity occurs at t = 2000, because the tank would be totally empty at this point. 28.6 A heat balance for the sphere can be written as dH = hA( Ta − T ) dt The heat gain can be transformed into a volume loss by considering the latent heat of fusion. Thus, dV hA =− (Ta − T ) dt ρL f (1) where ρ = density ≅ 1 kg/m3 and Lf = latent heat of fusion ≅ 333 kJ/kg. The volume and area of a sphere are computed by V= 4 3 πr 3 (2) A = 4πr 2 These can be combined with (1) to yield, dV = dt 3V h 4π   4 π  ρL f 2 /3 ( Ta − T ) This equation can be integrated along with the initial condition, V0 = 4 π ( 0.05) 3 = 0.000524 m 3 3 to yield the resulting volume as a function of time. V 0.0006 0.0003 t (s) 0 0 100000 200000 300000 This result can be converted into diameter using (2) d 0.1 0.05 t (s) 0 0 100000 200000 300000 28.7 The system for this problem is stiff. Thus, the use of a simple explicit Runge-Kutta scheme would involve using a very small time step in order to maintain a stable solution. A solver designed for stiff systems was used to generate the solution shown below. Two views of the solution are given. The first is for the entire solution domain. cb 400 200 ca 0 0 10 20 cc 30 In addition, we can enlarge the initial part of the solution to illustrate the fast transients that occur as the solution moves from its initial conditions to its dominant trajectories. cb 400 ca 200 0 0 cc 0.005 0.01 0.015 0.02 28.8 Several methods could be used to obtain a solution for this problem (e.g., finite-difference, shooting method, finite-element). The finite-difference approach is straightforward: D Ai −1 − 2 Ai + Ai +1 ∆x 2 − kAi = 0 Substituting parameter values and collecting terms gives − 1 × 10 −6 Ai −1 + (2 × 10 −6 −6 2 + 4 × 10 ∆ x ) − 1 × 10 −6 Ai +1 = 0 Using a ∆x = 0.2 cm this equation can be written for all the interior nodes. The resulting linear system can be solved with an approach like the Gauss-Seidel method. The following table and graph summarize the results. 0.1 0.05 0 0 1 2 3 4 28.9 The ODE to be solved is dP b A sin ωt =− P+ dt a a Substituting the parameters, it becomes dP = sin t − P dt The following Matlab script uses Euler’s method to solve the problem. dt=0.05; max=5; n=max/dt+1; t=zeros(1,n); p=zeros(1,n); t(1)=0; p(1)=90; for i=1:n p(i+1)=p(i)+dydt(t(i),p(i))*dt; t(i+1)=t(i)+dt; end plot(t,p) grid xlabel('Time-sec') ylabel('Pressure-mmHg') title('Pressure vs Time') zoom function s=dydt(t,p); A=1; w=1; s=A*sin(w*t)-p; P res s ure vs Tim e 90 80 70 P res s ure-m m Hg 60 50 40 30 20 10 0 -10 0 1 2 3 4 5 6 Tim e-s ec 28.10 Excel can be used to compute the basic results. As can be seen, the person died 1.13 hrs prior to being discovered. The non-self-starting Heun yielded the following time series of temperature: 120 100 80 60 0 1 2 3 28.11 The classical 4th order RK method yields 4 28.12 The classical 4th order RK method yields 28.13 (a) The first few steps of Euler’s method are shown in the following table A plot of the entire simulation is shown below: 8 4 0 0 5 10 15 20 Notice that because the Euler method is lower order, the peaks are increasing, rather than repeating in a stable manner as time progresses. This result is reinforced when a state-space plot of the calculation is displayed. 6 3 0 0 5 10 (b) The first few steps of the Heun method is shown in the following table A plot of the entire simulation is shown below: 8 4 0 0 5 10 15 20 Notice that in contrast to the Euler method, the peaks are stable manner as time progresses. This result is also reinforced when a state-space plot of the calculation is displayed. 4 2 0 0 5 (c) The first few steps of the 4th-order RK method is shown in the following table The results are quite close to those obtained with the Heun method in part (b). In fact, both the time series and state-space plots are indistinguishable from each other. 28.14 Using the step size of 0.1, (a) and (b) both give unstable results. The 4th-order RK method yields a stable solution. The first few values are shown in the following table. A plot of the result for x is also shown below. Notice how after about t = 6, this solution diverges from the double precision version in Fig. 28.9. 30 20 10 0 -10 0 5 10 15 20 -20 -30 28.15 The second-order equation can be reexpressed as a pair of first-order equations, dy =w dz dw f = ( L − z) 2 dz 2 EI We used Euler’s method with h = 1 to obtain the solution: z 30 20 10 y -0.5 0 0.5 1 28.16 The second-order equation can be reexpressed as a pair of first-order equations, dy =w dz dw 200 ze −2 z / 30 = ( L − z) 2 dz (5 + z) 2 EI We used Euler’s method with h = 1 to obtain the solution: z 30 20 10 y -0.5 0 0.5 1 28.17 This problem was solved using the Excel spreadsheet in a fashion similar to the last example in Sec. 28.1. We set up Euler’s method to solve the 3 ODEs using guesses for the diffusion coefficients. Then we formed a column containing the squared residuals between our predictions and the measured values. Adjusting the diffusion coefficients with the Solver tool minimized the sum of the squares. At first, we assumed that the diffusion coefficients were zero. For this case the Solver did not converge on a credible answer. We then made guesses of 1×107 for both. This magnitude was based on the fact that the volumes were of this order of magnitude. The resulting simulation did not fit the data very well, but was much better than when we had guessed zero. When we used Solver, it converged on E12 = 9.22×105 and E13 = 2.19×106 which corresponded to a sum of the squares of residuals of 2.007. Some of the Euler calculations are displayed below along with a plot of the fit. 100 50 0 0 5 10 15 20 It should be noted that we made up the “measurements” for this problem using the 4th-order RK method with values for diffusive mixing of E12 = 1×106 and E13 = 2×106. We then used a random number generator to add some error to this “data.” 28.18 The Heun method can be used to compute The results can be plotted. In addition, linear regression can be used to fit a straight line to lnp versus t to give lnp = 8.52 + 0.07t. Thus, as would be expected from a first-order model, the slope is equal to the growth rate of the population. p 20000 ln p 10 15000 9.6 10000 9.2 5000 8.8 0 8.4 0 10 20 0 10 20 28.19 The Heun method can be used to compute The results can be plotted. 20000 15000 10000 5000 0 0 10 20 The curve is s-shaped. This shape occurs because initially the population is increasing exponentially since p is much less than pmax. However, as p approaches pmax, the growth rate decreases and the population levels off. 28.20 (a) Nonlinear regression (e.g., using the Excel solver option) can be used to minimize the sum of the squares of the residuals between the data and the simulation. The resulting estimates are: a = 0.32823, b = 0.01231, c = 0.22445, and d = 0.00029. The fit is: 1200 moose 800 400 wolves 0 1960 (b) The results in state space are, 1990 2020 Wolves 300 200 100 0 0 500 1000 1500 Moose (c) 1200 800 400 0 1960 1990 2020 Wolves 300 200 100 0 0 500 1000 1500 Moose (d) 1200 800 400 0 1960 1990 2020 Wolves 300 200 100 0 0 500 1000 1500 Moose 28.21 Main Program: % Hanging static cable - w=w(x) % Parabolic solution w=w(x) % CUS Units (lb,ft,s) % w = wo(1+sin(pi/2*x/l) es=0.5e-7 % Independent Variable x, xs=start x, xf=end x xs=0; xf=200; %initial conditions [y(1)=cable y-coordinate, y(2)=cable slope]; ic=[0 0]; global wToP wToP=0.0025; [x,y]=ode45('slp',xs,xf,ic,.5e-7); yf(1)=y(length(x)); wTo(1)=wToP; ea(1)=1; wToP=0.002; [x,y]=ode45('slp',xs,xf,ic,.5e-7); yf(2)=y(length(x)); wTo(2)=wToP; ea(2)=abs( (yf(2)-yf(1))/yf(2) ); for k=3:10 wTo(k)=wTo(k-1)+(wTo(k-1)-wTo(k-2))/(yf(k-1)-yf(k-2))*(50-yf(k-1)); wToP=wTo(k); [x,y]=ode45('slp',xs,xf,ic,.5e-7); yf(k)=y(length(x)); ea(k)=abs( (yf(k)-yf(k-1))/yf(k) ); if (ea(k)<=es) %Analytic Solution with constant w (for Comparison) xa=xs:.01:xf; ya=(0.00125)*(xa.*xa); plot(x,y(:,1),xa,ya,'--'); grid; xlabel('x-coordinate - ft'); ylabel('y-coordinate - ft'); title('Cable - w=wo(1+sin(pi/2*x/l))'); fprintf('wTo %f\n', wTo) fprintf('yf %f\n', yf) fprintf('ea %f\n', ea) break end end Function ‘slp’: function dxy=slp(x,y) global wToP dxy=[y(2);(wToP)*(1+sin(pi/2*x/200))] Cable - w= wo(1+ s in(pi/2*x /l)) 50 45 40 y -c oordinate - ft 35 30 25 20 15 10 5 0 0 20 40 60 80 100 120 x -c oordinate - ft 28.22 140 160 180 200 28.23 The second-order equation can be reexpressed as a pair of first-order equations, dq =i dt di = −0.05i − 4q + sin 18708 . t dz The parameters can be substituted and the system solved with the 4th-order RK method in double-precision with h = 0.1. A table showing the first few steps and a graph of the entire solution are shown below. 4 2 0 -2 -4 0 20 40 60 80 100 28.24 The second-order equation can be reexpressed as a pair of first-order equations, dq =i dt di R q =− i− dt L CL The parameters can be substituted and the system solved with the 4th-order RK method with h = 0.005. A table showing the first few steps and a graph of the entire solution are shown below. i 4 2 0 -2 0 0.1 0.2 0.3 0.4 0.5 -4 -6 q 0.08 0.04 0 -0.04 0 0.1 0.2 -0.08 28.25 The equation can be solved analytically as 0.3 0.4 0.5 di R =− i dt L di R = − dt i L ln i = −( R / L ) t + C C = 0.001 i = 0.001e −2t The numerical solution can be obtained by expressing the equation as di = −2i dt and using Euler’s method with h = 0.05 to solve for the current. Some selected values are shown, along with a plot of both the analytical and numerical result. A better match would be obtained by using a smaller step or a higher-order method. 0.6 analytical 0.4 Euler 0.2 0.0 0 0.5 1 28.26 The numerical solution can be obtained by expressing the equation as di = −( −i + i 3 ) 2 dt and using Euler’s method with h = 0.05 to solve for the current. Some selected values are shown, along with a plot of the numerical result. Note that the table and plot also show the analytical solution for the linear case computed in Prob. 28.19. 0.6 analytical 0.4 Euler 0.2 0.0 0 0.5 1 28.27 Using an approach similar to Sec. 28.3, the system can be expressed in matrix form as 1 − λ  − 1 − 1  i1  = {0} 2 − λ  i 2  A package like MATLAB can be used to evaluate the eigenvalues and eigenvectors as in >> a=[1 -1;-1 2]; >> [v,d]=eig(a) v = d = 0.8507 0.5257 -0.5257 0.8507 0.3820 0 0 2.6180 Thus, we can see that the eigenvalues are λ = 0.382 and 2.618 or natural frequencies of ω = 0.618/ LC and 1.618/ LC . The eigenvectors tell us that these correspond to oscillations that coincide (0.8507 0.5257) and which run counter to each other (−0.5257 0.8507). 28.28 The differential equations to be solved are linear: nonlinear: dθ =v dt dv 32.2 =− θ dt 4 dθ =v dt dv 32.2 =− sin θ dt 4 A few steps for the 4th-order RK solution of the nonlinear system are contained in the following table and a plot of both solutions is shown below. lin-theta 6 lin-v nonlin-theta nonlin-v 3 0 0 2 4 6 -3 28.29 The differential equations to be solved are dx =v dt dv c k =− v− v dt m m A few steps for the 2nd-order RK solution (Heun without iteration) are shown in the following table and a plot of displacement is shown below. 0.2 0 0 0.1 0.2 0.3 0.4 -0.2 28.30 The differential equation to be solved is dT = 0.2 (40 − T ) dt A few steps for the 2nd-order RK solution (Heun without iteration) are shown in the following table and a plot of temperature versus time is shown below. The temperature will drop 95% of the way to the new temperature in 3/0.2 = 15 minutes. 100 50 0 0 5 10 15 20 25 28.31 The differential equation to be solved is dQ 100 (20 − 2.5)(20 − t ) = 0.4(10) dt 100 − 2.5t A few steps for the 2nd-order RK solution (Heun without iteration) are shown in the following table and a plot of heat flow versus time is shown below. 20000 10000 0 0 5 10 15 20 25 28.32 The differential equations to be solved are nonlinear: dv 0.235 2 = 9.8 − v dt 68.1 linear: dv 12.5 = 9.8 − v dt 68.1 A few steps for the solution (Euler) are shown in the following table, which also includes the analytical solution from Example 1.1. A plot of the result is also shown below. Note, the nonlinear solution is the bolder line. 60 40 20 0 0 10 20 30 28.33 The differential equations to be solved are dv 12.5 = 9.8 − v dt 68.1 t < 15 s: t ≥ 15 s: dv 50 = 9.8 − v dt 68.1 The first few steps for the solution (Euler) are shown in the following table, along with the steps when the parachute opens. A plot of the result is also shown below. 60 40 20 0 0 28.34 10 20 30 %Damped spring mass system %mass: m=1 kg %damping, nonlinear: a sgn(dx/dt) (dx/dt)^2, a=2 N/(m/s)^2 %spring, nonlinear: bx^3, b=5 N/m^3 % MATLAB 5 version %Independent Variable t, tspan=[tstart tstop] %initial conditions [x(1)=velocity, x(2)=displacement]; t0=0; tf=8; tspan=[0 8]; ic=[1 0.5]; % a) linear solution [t,x]=ode45('kc',tspan,ic); subplot(221) plot(t,x); grid; xlabel('time - sec.'); ylabel('displacement - m; velocity - m/s'); title('d2x/dt2+2(dx/dt)+5x=0') subplot(222) %Phase-plane portrait plot(x(:,2),x(:,1)); grid; xlabel('displacement - m'); ylabel('velocity - m/s'); title('d2x/dt2+2(dx/dt)+5x=0'); % b) nonlinear spring [t,x]=ode45('nlk',tspan,ic); subplot(223) plot(t,x); grid; xlabel('time - sec.'); ylabel('displacement - m; velocity - m/s'); title('d2x/dt2+2(dx/dt)+5x^3=0') %Phase-plane portrait subplot(224) plot(x(:,2),x(:,1)); grid; xlabel('displacement - m'); ylabel('velocity - m/s'); title('d2x/dt2+2(dx/dt)+5x=0'); pause % c) nonlinear damping [t,x]=ode45('nlc',tspan,ic); subplot(221) plot(t,x); grid; xlabel('time - sec.'); ylabel('displacement - m; velocity - m/s'); title('d2x/dt2+2sign(dx/dt)(dx/dt)^2+5x=0') %Phase-plane portrait subplot(222) plot(x(:,2),x(:,1)); grid; xlabel('displacement - m'); ylabel('velocity - m/s'); title('d2x/dt2+2sign(dx/dt)(dx/dt)^2+5x=0'); % d) nonlinear damping and spring [t,x]=ode45('nlck',tspan,ic); subplot(223) plot(t,x); grid; xlabel('time - sec.'); ylabel('displacement - m; velocity - m/s'); title('d2x/dt2+2sign(dx/dt)(dx/dt)^2+5x^3=0') %Phase-plane portrait subplot(224) plot(x(:,2),x(:,1)); grid; xlabel('displacement - m'); ylabel('velocity - m/s'); title('d2x/dt2+2sign(dx/dt)(dx/dt)^2+5x^3=0'); Functions: %Damped spring mass system - m d2x/dt2 + c dx/dt + k x =0 %mass: m=1 kg % linearc=2 N/(m/s) % lineark=5 N/m %x(1)=velocity, x(2)=displacement function dx=kc(t,x); dx=[-2*x(1)-5*x(2); x(1)] %Damped spring mass system - m d2x/dt2 + c dx/dt + k x =0 %mass: m=1 kg %damping: linearc=2 N/(m/s) %spring: nonlinear- kx=bx^3, b=5 N/m^3 function dx=nlk(t,x); dx=[-2*x(1)-5*x(2).*x(2).*x(2); x(1)] %Damped spring mass system - m d2x/dt2 + c dx/dt + k x =0 %mass: m=1 kg %damping: nonlinear- c dx/dt = a sgn(dx/dt) (dx/dt)^2, a=2 N/(m/s)^2 %spring: linearkx=5x %x(1)=velocity, x(2)=dispacement function dx=nlc(t,x); dx(1)=-2*sign(x(1))*x(1)*x(1)-5*x(2); dx(2)= x(1); %Damped spring mass system - m d2x/dt2 + c dx/dt + k x =0 %mass: m=1 kg %damping: nonlinear- c dx/dt = a sgn(dx/dt) (dx/dt)^2, a=2 N/(m/s)^2 %spring: nonlinear- k x = bx^3, b=5 N/m^3 %x(1)=velocity, x(2)=dispacement function dx=nlck(t,x); dx=[-2*sign(x(1)).*x(1).*x(1)-5*x(2).*x(2).*x(2); x(1)] d2x /dt2+ 2(dx /dt)+ 5x = 0 0.5 0.5 veloc ity - m /s 1 0 -0.5 -1 0 2 4 6 0 -0.5 -1 -0.5 8 tim e - s ec . d2x /dt2+ 2(dx /dt)+ 5x 3 = 0 0 0.5 1 dis plac em ent - m d2x /dt2+ 2(dx /dt)+ 5x = 0 1 veloc ity - m /s 1 0.5 0 -0.5 0.5 0 -0.5 0 2 4 6 8 0 tim e - s ec . d2x /dt2+ 2s ign(dx /dt)(dx /dt) 2 + 5x = 0 1 0 -0.5 -1 2 4 6 0.5 0 -0.5 -1 2 4 tim e - s ec . 0.6 0.8 6 8 0.5 0 -0.5 -1 -0.5 8 tim e - s ec . d2x /dt2+ 2s ign(dx /dt)(dx /dt) 2 + 5x 3 = 0 1 0 0.4 d2x /dt2+ 2s ign(dx /dt)(dx /dt) 2 + 5x = 0 1 veloc ity - m /s 0.5 0 0.2 dis plac em ent - m 0 0.5 1 dis plac em ent - m d2x /dt2+ 2s ign(dx /dt)(dx /dt) 2 + 5x 3 = 0 1 veloc ity - m /s dis plac em ent - m ; veloc ity - m /s dis plac em ent - m ; veloc ity - m /s dis plac em ent - m ; veloc ity - m /s dis plac em ent - m ; veloc ity - m /s 28.35 d2x /dt2+ 2(dx /dt)+ 5x = 0 1 0.5 0 -0.5 -1 -0.5 0 0.5 dis plac em ent - m %Forced damped spring-mass system w/ material damping %mass: m=2 kg %damping, nonlinear material: b sgn(dx/dt) abs(x), b=1 N/m %spring, linear: kx = 6x %forcing function: F=Fo(sin(wt)), Fo=2 N, w=0.5 rad/s 1 % MATLAB 5 version %Independent Variable t, tspan=[tstart tstop] %initial conditions [x(1)=velocity, x(2)=displacement]; tspan=[0 15]; ic=[0 1]; [t,x]=ode45('nlF',tspan,ic); ts=0:.01:15; Sin=2*sin(0.5*ts); plot(t,x,ts,Sin,'--'); grid; xlabel('time - sec.'); ylabel('displacement - m; velocity - m/s; force - N'); title('non-linear, forced,damped spring-mass system, time response') Function ‘n1F’: %Forced damped spring-mass system w/ material damping %mass: m=2 kg %damping, nonlinear air: b sgn(dx/dt) (dx/dt)^2, b=1 N/m %spring, linear: kx = 6x %forcing function: F=Fo(sin(wt), Fo=2 N, w=0.5 rad/s % x(1)= velocity, x(2)= displacement function dx=nlF(t,x); dx=[-0.5*sign(x(1)).*x(1).*x(1)-3*x(2)+sin(0.5*t); x(1)] non-linear, forc ed,dam ped s pring-m as s s y s tem , tim e res pons e 2 dis plac em ent - m ; veloc ity - m /s ; forc e - N 1.5 1 0.5 0 -0.5 -1 -1.5 -2 0 5 10 15 tim e - s ec . 28.36 % % % % % % % ODE Boundary Value Problem Tapered conical cooling fin u[xx]+(2/x)(u[x]-pu)=0 BC. u(x=0)=0 u(x=1)=1 i=spatial index, from 1 to R numbering for points is i=1 to i=R for (R-1) dx spaces u(i=1)=0 and u(i=R)=1 R=21; %Constants dx=1/(R-1); dx2=dx*dx; %Parameters p(1)=10; p(2)=20; p(3)=50; p(4)=100; %sizing matrices u=zeros(1,R); x=zeros(1,R); a=zeros(1,R); b=zeros(1,R); c=zeros(1,R); d=zeros(1,R); ba=zeros(1,R); ga=zeros(1,R); %Independent Variable x=0:dx:1; %Boundary Conditions u(1)=0; u(R)=1; for k=1:4; %/Coefficients b(2)=-2-2*p(k)*dx2/dx; c(2)=2; for i=3:R-2, a(i)=1-dx/(dx*(i-1)); b(i)=-2-2*p(k)*dx2/(dx*(i-1)); c(i)=1+1/(i-1); end a(R-1)=1-dx/(dx*(R-2)); b(R-1)=-2-2*p(k)*dx2/(dx*(R-2)); d(R-1)=-(1+1/(R-2)); %Solution by Thomas Algorithm ba(2)=b(2); ga(2)=d(2)/b(2); for i=3:R-1, ba(i)=b(i)-a(i)*c(i-1)/ba(i-1); ga(i)=(d(i)-a(i)*ga(i-1))/ba(i); end %back substitution step u(R-1)=ga(R-1); for i=R-2:-1:2, u(i)=ga(i)-c(i)*u(i+1)/ba(i); end %Plot plot(x,u) title('u[xx]+(2/x)(u[x]-pu)=0; u(x=0)=0, u(x=1)=1') xlabel('x -ND Length') ylabel('u - ND Temperature') hold on end grid hold off gtext('p=10');gtext('p=20');gtext('p=50');gtext('p=100'); u[x x ]+ (2/x )(u[x ]-pu)= 0; u(x = 0)= 0, u(x = 1)= 1 1 0.9 0.8 u - ND Tem perature 0.7 0.6 0.5 0.4 p= 10 0.3 0.2 p= 20 0.1 p= 50 p= 100 0 0 0.1 0.2 0.3 0.4 0.5 0.6 x -N D Length t 0 0.1 0.2 0.3 0.4 0.5 v 0 0.98 1.942012 2.886365 3.813385 4.723389 dvdt 9.8 9.620117 9.443537 9.270197 9.100039 8.933005 0.7 0.8 0.9 1 • • • 14.9 15 15.1 15.2 15.3 15.4 15.5 x 0.3 0.299844 0.299379 0.298609 0.297539 50.01245 50.07445 47.37791 44.87936 42.56425 40.41912 38.43149 0.620036 -26.9654 -24.9855 -23.1511 -21.4513 -19.8763 -18.417 t 0 0.1 0.2 0.3 0.4 0.5 linear v 0 0.98 1.942012 2.886365 3.813385 4.723389 dvdt 9.8 9.620117 9.443537 9.270197 9.100039 8.933005 nonlinear v dvdt 0 9.8 0.98 9.796686 1.959669 9.786748 2.938343 9.770206 3.915364 9.747099 4.890074 9.717481 t 0 0.1 0.2 0.3 0.4 0.5 T 0 160.0003 319.6 478.7972 637.5899 795.976 k11 1598 1593.995 1589.97 1585.924 1581.859 1577.772 T-end 159.8 319.3997 478.597 637.3897 795.7758 953.7532 k21 1602.005 1598 1593.975 1589.929 1585.863 1581.777 phi1 1600.003 1595.997 1591.972 1587.927 1583.861 1579.774 t 0 0.1 0.2 0.3 0.4 0.5 T 90 89.01 88.0396 87.08842 86.15607 85.24218 k11 -10 -9.802 -9.60792 -9.41768 -9.23121 -9.04844 T-end 89 88.0298 87.07881 86.14665 85.23295 84.33733 k21 -9.8 -9.60596 -9.41576 -9.22933 -9.04659 -8.86747 phi1 -9.9 -9.70398 -9.51184 -9.32351 -9.1389 -8.95795 v 0 -0.31068 -0.61743 -0.91998 -1.21806 k11 0 -0.31068 -0.61743 -0.91998 -1.21806 k12 -312.5 -308.713 -304.65 -300.318 -295.726 x 0.3 0.299533 0.298761 0.297689 0.296321 analytical 0 0.971061 1.92446 2.860518 3.779552 4.681871 v -0.3125 -0.61939 -0.92208 -1.2203 -1.51379 k21 -0.3125 -0.61939 -0.92208 -1.2203 -1.51379 k22 -308.854 -304.787 -300.452 -295.856 -291.007 phi1 -0.15625 -0.46503 -0.76975 -1.07014 -1.36593 phi2 -310.677 -306.75 -302.551 -298.087 -293.366 t thet v k11 k12 thet v k21 k22 thet v k31 k32 thet v k41 k42 phi1 phi2 0 0.785 0.000 0.000 -5.692 0.785 -0.028 -0.028 -5.692 0.785 -0.028 -0.028 -5.691 0.785 -0.057 -0.057 -5.691 -0.028 -5.692 0.01 0.785 -0.057 -0.057 -5.691 0.785 -0.085 -0.085 -5.689 0.785 -0.085 -0.085 -5.688 0.784 -0.114 -0.114 -5.686 -0.085 -5.688 0.02 0.784 -0.114 -0.114 -5.686 0.784 -0.142 -0.142 -5.682 0.784 -0.142 -0.142 -5.682 0.783 -0.171 -0.171 -5.678 -0.142 -5.682 0.03 0.783 -0.171 -0.171 -5.678 0.782 -0.199 -0.199 -5.673 0.782 -0.199 -0.199 -5.672 0.781 -0.227 -0.227 -5.666 -0.199 -5.672 0.04 0.781 -0.227 -0.227 -5.666 0.780 -0.256 -0.256 -5.660 0.780 -0.256 -0.256 -5.659 0.778 -0.284 -0.284 -5.652 -0.256 -5.659 t (analytical) 0 0.600000 0.05 0.542902 0.1 0.491238 0.15 0.444491 0.2 0.402192 0.25 0.363918 (Euler) 0.600000 0.561600 0.523153 0.485155 0.448059 0.412248 didt -0.768000 -0.768949 -0.759943 -0.741923 -0.716216 -0.684375 t (analytical) 0 0.600000 0.05 0.542902 0.1 0.491238 0.15 0.444491 0.2 0.402192 0.25 0.363918 t 0 0.005 0.01 0.015 i -3.282 -3.837 -4.157 -4.250 (Euler) 0.600000 0.540000 0.486000 0.437400 0.393660 0.354294 didt -1.200000 -1.080000 -0.972000 -0.874800 -0.787320 -0.708588 q k11 k12 imid qmid k21 k22 imid qmid k31 k32 iend 0.100 -134.37 -3.282 -3.617 0.092 -111.24 -3.617 -3.560 0.091 -110.7 -3.560 -3.835 0.082 -87.485 -3.837 -4.055 0.073 -63.93 -4.055 -3.996 0.072 -64.01 -3.996 -4.157 0.062 -40.917 -4.157 -4.259 0.052 -18.09 -4.259 -4.202 0.051 -18.72 -4.202 -4.251 0.041 3.159 -4.250 -4.242 0.030 24.25 -4.242 -4.189 0.030 23.16 -4.189 -4.134 0.02 -4.133 0.020 t i 0 0.000 0 0.1 0.009 3 0.2 0.036 4 0.3 0.079 1 0.4 0.134 2 42.892 -4.133 -4.025 0.010 q 0.000 0 0.000 3 0.002 5 0.008 1 0.018 7 k11 0.000 0 0.184 3 0.353 9 0.495 8 0.598 9 k12 0.000 0 0.009 3 0.036 4 0.079 1 0.134 2 imid 0.000 0 0.018 5 0.054 1 0.103 9 0.164 2 61.41 -4.025 -3.979 0.010 qmid 0.000 0 0.000 8 0.004 3 0.012 1 0.025 4 k21 0.093 4 0.272 9 0.431 1 0.555 5 0.636 1 k22 0.000 0 0.018 5 0.054 1 0.103 9 0.164 2 imid 0.004 7 0.022 9 0.057 9 0.106 9 0.166 0 qmid 0.000 0 0.001 2 0.005 2 0.013 3 0.026 9 59.95 -3.979 -3.833 k31 0.093 2 0.270 9 0.427 3 0.550 4 0.630 0 k32 0.004 7 0.022 9 0.057 9 0.106 9 0.166 0 t 0 0.5 1 1.5 • • • 18 18.5 19 19.5 20 p 5000 5384.023 5786.181 6205.641 k1 750 786.9276 822.4373 856.0284 pend 5375 5777.487 6197.4 6633.655 k2 786.0938 821.7039 855.4023 886.6772 phi 768.0469 804.3157 838.9198 871.3528 18480.96 18615.33 18738.61 18851.6 18955.02 280.733 257.7616 236.3663 216.4921 198.0754 18621.33 18744.21 18856.8 18959.84 19054.06 256.7271 235.3885 215.5715 197.2119 180.2397 268.7301 246.575 225.9689 206.852 189.1575 t 0 0.5 1 1.5 2 • • • 18 18.5 19 19.5 20 p 5000 5178.063 5362.466 5553.437 5751.209 k1 350 362.4644 375.3726 388.7406 402.5846 pend 5175 5359.295 5550.153 5747.807 5952.501 k2 362.25 375.1506 388.5107 402.3465 416.6751 phi 356.125 368.8075 381.9417 395.5436 409.6299 17622.69 18250.28 18900.22 19573.3 20270.36 1233.588 1277.52 1323.015 1370.131 1418.925 18239.48 18889.04 19561.72 20258.37 20979.82 1276.764 1322.233 1369.321 1418.086 1468.587 1255.176 1299.876 1346.168 1394.108 1443.756 t 0 0.1 c1 0 2.193515 c2 0 0 c3 100 95.61297 dc1/dt 21.93515 19.41211 dc2/dt 0 0.252723 iend 0.009 3 0.036 4 0.079 1 0.134 2 0.197 2 qend k41 k42 phi1 phi2 0.082 -87.70 -3.835 -111.0 -3.578 0.062 -41.12 -4.157 -64.08 -4.016 0.041 2.976 -4.251 -18.59 -4.222 0.020 42.736 -4.134 23.45 -4.208 2 0.000 76.688 -3.833 60.38 -3.996 3 qend 0.000 5 0.002 6 0.008 2 0.018 8 0.035 3 dc3/dt -43.8703 -40.9834 k41 0.183 7 0.353 3 0.495 3 0.598 5 0.653 8 k42 0.009 3 0.036 4 0.079 1 0.134 2 0.197 2 phi1 0.092 8 0.270 9 0.427 7 0.551 0 0.630 8 phi2 0.0031 0.0214 0.0566 0.1058 0.1653 0.2 0.3 0.4 0.5 4.134726 5.848151 7.356012 8.678451 0.025272 0.072619 0.139161 0.222309 91.51463 87.68125 84.09121 80.72481 17.13425 15.07861 13.22439 11.55268 0.473465 0.66542 0.83148 0.974263 z 0 1 2 3 4 5 • • • 26 27 28 29 30 y 0 0 0.0036 0.010564 0.020664 0.03368 w 0 0.0036 0.006964 0.0101 0.013016 0.01572 dydz 0 0.0036 0.006964 0.0101 0.013016 0.01572 dwdz 0.0036 0.003364 0.003136 0.002916 0.002704 0.0025 0.676 0.7137 0.751464 0.789264 0.82708 0.0377 0.037764 0.0378 0.037816 0.03782 0.0377 0.037764 0.0378 0.037816 0.03782 0.000064 0.000036 0.000016 0.000004 0 z 0 1 2 3 4 5 • • • 26 27 28 29 30 f(z) 0 31.18357 50.0099 61.40481 68.08252 71.65313 y 0 0 0 0.002098 0.007333 0.016148 w 0 0 0.002098 0.005235 0.008816 0.012498 dydz 0 0 0.002098 0.005235 0.008816 0.012498 dwdz 0 0.002098 0.003137 0.003581 0.003682 0.003583 29.63907 27.89419 26.24164 24.67818 23.20033 0.700979 0.741693 0.782444 0.823216 0.863995 0.040713 0.040751 0.040771 0.04078 0.040782 0.040713 0.040751 0.040771 0.04078 0.040782 3.79E-05 2.01E-05 8.4E-06 1.97E-06 0 z 0 1 2 3 4 5 • • • 26 27 28 29 30 y 0 0 0.0036 0.010564 0.020664 0.03368 w 0 0.0036 0.006964 0.0101 0.013016 0.01572 dydz 0 0.0036 0.006964 0.0101 0.013016 0.01572 dwdz 0.0036 0.003364 0.003136 0.002916 0.002704 0.0025 0.676 0.7137 0.751464 0.789264 0.82708 0.0377 0.037764 0.0378 0.037816 0.03782 0.0377 0.037764 0.0378 0.037816 0.03782 0.000064 0.000036 0.000016 0.000004 0 t 0 0.1 0.2 x 5 9.78147 17.70297 y 5 17.07946 20.8741 Z 5 10.43947 35.89688 -38.3338 -35.9004 -33.664 -31.607 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10.81088 0.549577 -3.16461 -5.57588 -8.88719 -11.9142 -10.6668 -6.84678 -2.52924 -5.54419 -5.84129 -8.42037 -12.6789 -13.43 -7.21784 -3.43018 0 0.1 0.2 0.3 0.4 0.5 2 1.887095 1.787897 1.701588 1.627287 1.564109 3 2.935517 2.863301 2.785107 2.702536 2.617016 0 0.1 0.2 0.3 0.4 0.5 2 1.886984 1.787729 1.701406 1.627125 1.56399 3 2.935308 2.862899 2.784535 2.701821 2.616185 t 0 0.1 0.2 0.3 0.4 0.5 x 2 1.88 1.773968 1.681301 1.601231 1.532907 y 3 2.94 2.870616 2.793738 2.711153 2.624496 t 0 0.0625 0.125 0.1875 0.25 0.3125 0.375 0.4375 0.5 0.5625 0.625 0.6875 0.75 0.8125 0.875 0.9375 1 • • • x 0 0.2 0.4 C 1 0.941218 0.885749 0.833497 0.784309 0.738024 0.694475 0.653506 0.614963 0.578703 0.54459 0.512497 0.482304 0.453896 0.427168 0.40202 0.378358 A 0.1 0.067208 0.045169 39.30744 28.07461 22.36888 19.92312 22.14149 29.80001 33.39903 29.30716 Te 25 66.18648 85.80247 93.93385 96.02265 94.99472 92.41801 89.12894 85.57041 81.97385 78.45733 75.07829 71.86194 68.81648 65.9413 63.23134 60.67946 x A x A x A 1.2 1.4 0.009215 0.006193 2.2 2.4 0.001263 0.000848 3.2 3.4 0.000166 0.000106 0.6 0.8 1 0.030357 0.020402 0.013712 t 0 5 10 15 20 25 M 1000 997.5 995 992.5 990 987.5 t 0 2 4 6 8 10 t 0 1 2 3 4 5 • • • 76 77 78 79 80 t 0 10 20 30 40 50 1.6 1.8 2 cin 0 9.063462 16.484 22.55942 27.53355 31.60603 m 8000 8000 7998.997 7997.038 7994.164 7990.419 c 10 9.503173 9.832427 10.7681 12.13698 13.80328 c1 0.0000 1.0000 1.8600 2.5996 3.2357 3.7827 c2 c3 0.0000 0.0000 0.0000 5.0000 0.2000 8.6250 0.5320 11.2581 0.9455 13.1754 1.4035 14.5758 7.1428 7.1428 7.1428 7.1428 7.1428 7.1426 7.1426 7.1427 7.1427 7.1427 c 10 25.72917 35.27317 41.06419 44.57801 46.71009 k1 2 1.213542 0.736342 0.446791 0.2711 0.164495 18.8311 18.8311 18.8311 18.8311 18.8311 0.004162 0.002797 0.00188 s 8 8.02005 8.039193 8.057469 8.074914 8.091563 k1 -0.5 -0.02199 0.332579 0.589566 0.769829 0.890138 2.6 2.8 3 0.000569 0.00038 0.000253 3.6 3.8 4 6.23E-05 2.88E-05 0 dmdt 0 -0.2005 -0.39193 -0.57469 -0.74914 -0.91563 cend 9 9.459202 10.49758 11.94723 13.67664 15.58355 cin-end 9.063462 16.484 22.55942 27.53355 31.60603 34.94029 k2 0.003173 0.35124 0.603092 0.779316 0.89647 0.967837 phi -0.24841 0.164627 0.467835 0.684441 0.833149 0.928987 c4 0.0000 0.0000 0.5625 1.4351 2.4528 3.5102 c5 0.0000 0.0000 0.0300 0.0900 0.1785 0.2933 dc1dt 1.0000 0.8600 0.7396 0.6361 0.5470 0.4704 dc2dt 0.0000 0.2000 0.3320 0.4135 0.4580 0.4758 dc3dt 5.0000 3.6250 2.6331 1.9173 1.4004 1.0267 dc4dt 0.0000 0.5625 0.8726 1.0176 1.0575 1.0328 dc5dt 0.0000 0.0300 0.0600 0.0885 0.1147 0.1380 13.0962 13.0980 13.0997 13.1013 13.1028 7.0053 7.0135 7.0213 7.0286 7.0354 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0018 0.0017 0.0016 0.0015 0.0014 0.0082 0.0078 0.0073 0.0069 0.0064 cmid 20 31.79688 38.95487 43.29814 45.93351 47.53257 k2 1.5 0.910156 0.552256 0.335093 0.203325 0.123371 t c k1 c 0 10 2 30 10 25 1.25 37.5 20 34.375 0.78125 42.1875 30 40.23438 0.488281 45.11719 40 43.89648 0.305176 46.94824 50 46.1853 0.190735 48.09265 cmid 17.5 30.27995 38.03445 42.73965 45.59463 47.32695 k2 1 0.625 0.390625 0.244141 0.152588 0.095367 k3 1.625 0.986003 0.598278 0.363017 0.220268 0.133652 phi 1.5 0.9375 0.585938 0.366211 0.228882 0.143051 cend 26.25 35.58919 41.25594 44.69436 46.78069 48.04662 k4 1.1875 0.72054 0.437203 0.265282 0.160965 0.097669 phi 1.572917 0.9544 0.579102 0.351382 0.213208 0.129369