Kartho Matlab.pdf

1

To Visualize 2D and 3D Plots Aim: To visualize 2D and 3D plots using MATLAB with the help of Aerospace related questions and interpreting it’s answers.

Theory: A two-dimensional graph is a set of points in two-dimensional space. If the points are real and if Cartesian coordinates are used, each axis depicts the potential values of a particular real variable. 3D plot are generated from data defined as Z=f(X,Y). As for 2D plots, there are two ways to obtain a 3D plot depending on the way the (X,Y,Z) values are defined. Types of plot: 1. Line Plot: Line plots are a useful way to compare sets of data or track changes over time. 2D or 3D line plot can be drafted on linear or logarithmic scale based on the requirements.

2. Data Distribution Plots: Visualization of the distribution of data using plots such as histograms, pie charts, or word clouds is known as Data Distribution Plots. 3. Discrete Data Plots: Visualizing discrete data using plots such as bar graphs or stem plots. For example, one can create a vertical or horizontal bar graph where the bar lengths are proportional to the values that they represent.

4. Geographical Plots: It is a plot of the movements of one or more craft relative to the surface of the Earth.

5. Polar Plots: The are the data plots which are done in polar Coordinate System.

6. Contour Plots: A contour plot is a geographical technique for representing a 3-D surface by plotting constant z slices, called the contours, on a 2 Dimensional format.

7. Vector Fields: a vector field is an assignment of a vector to each point in a subset of space. A vector field in the plane can be visualised as a collection of arrows with a given magnitude and direction, each attached to a point in the plane.

2

8. Animation Plots: Animations are generated from a list (or other iterable) of graphics objects.

Questions: 1. Consider a gulfstream 4 twin-turbofan executive transport aircraft. Calculate and plot the thrust required curve at an altitude of 13000 ft. Assuming weight of the aircraft is 73000 pounds. The airplane data are S=950 sq ft, AR= 5.92, Cdo=0.015, k=0.08.

First, we select a range of Velocities required for the graph between thrust and velocity. Next we calculate coefficient of lift using the formula: Cl= 2(WEIGHT)/(DESITY)*(VELOCITY)^2*(AREA)

Next we calculate coefficient of drag using the formula: Cd=Cdo+k*Cl^2 Which gives: Cd= 0.015+0.08*Cl^2

Next Thrust is calculated with the formula: Thrust=0.5*(DENSITY)*(VELOCITY)^2*(AREA)*(Cd) And thus graph between thrust and velocity is plotted.

MATLAB Results: We put the following code in matlab: velocity=100:100:1300 density=0.00089068 weight=73000 area=950 Coeff_lift=2*weight./density*velocity.^2*area Coeff_drag=0.015+0.08*Coeff_lift.^2 Thrust=0.5.*density.*velocity.^2.*area.*Coeff_drag drag=0.5.*density.*velocity.^2.*area.*Coeff_lift

3

In command window: plot(velocity,Thrust)

Result in command window: >> velocity =

Columns 1 through 11

100 1000

200

300

1100

Columns 12 through 13

1200

density =

8.9068e-04

weight =

73000

area =

1300

400

500

600

700

800

900

4 950

Coeff_lift =

1.0e+17 *

0.0156 0.0623 0.1402 0.2492 1.5572 1.8843 2.2424 2.6317

0.3893

0.5606

0.7630

0.9966

1.2614

0.1212

0.2514

0.4658

0.7946

1.2728

0.0128

0.0383

0.0966

0.2152

0.4362

Coeff_drag =

1.0e+33 *

0.0002 0.0031 0.0157 0.0497 1.9400 2.8403 4.0228 5.5408

Thrust =

1.0e+39 *

0.0000 0.0001 0.0006 0.0034 0.8208 1.4540 2.4508 3.9616

drag =

5 1.0e+23 *

0.0001 0.0011 0.0053 0.0169 0.6588 0.9646 1.3661 1.8817

Graph Obtained:

0.0412

0.0854

0.1582

0.2699

0.4323

6

2. From the above Values, plot and analyse variation of Cl^(3/2)/Cd, Cl/Cd, Cl^(1/2)/Cd Curves.

velocity=100:100:1300 density=0.00089068 weight=73000 area=950 Coeff_lift=2*weight./density*velocity.^2*area Coeff_drag=0.015+0.08*Coeff_lift.^2 R1=Coeff_lift.^(3/2) Ratio1=R1./Coeff_drag Ratio2=Coeff_lift./Coeff_drag R2=Coeff_lift.^0.5 Ratio3=R2./Coeff_drag

Program Compilation: >> KarthoVampee2

velocity =

Columns 1 through 5

100

200

300

400

500

800

900

1000

Columns 6 through 10

600

700

Columns 11 through 13

1100

density =

1200

1300

7

8.9068e-04

weight =

73000

area =

950

Coeff_lift =

1.0e+17 *

Columns 1 through 6

0.0156

0.0623

0.1402

0.2492

0.3893

0.5606

1.5572

1.8843

2.2424

Columns 7 through 12

0.7630

0.9966

Column 13

2.6317

1.2614

8

Coeff_drag =

1.0e+33 *

Columns 1 through 6

0.0002

0.0031

0.0157

0.0497

0.1212

0.2514

1.9400

2.8403

4.0228

0.0393

0.0768

0.1327

0.6145

0.8179

1.0619

Columns 7 through 12

0.4658

0.7946

1.2728

Column 13

5.5408

R1 =

1.0e+26 *

Columns 1 through 6

0.0006

0.0049

0.0166

Columns 7 through 12

0.2108

0.3146

0.4480

9

Column 13

1.3501

Ratio1 =

1.0e-06 *

Columns 1 through 6

0.3168 0.1584

0.1056

0.0792

0.0634

0.0528

0.0317

0.0288

0.0264

0.0502

0.0321

0.0223

Columns 7 through 12

0.0453

0.0396

0.0352

Column 13

0.0244

Ratio2 =

1.0e-14 *

Columns 1 through 6

0.8027

0.2007

0.0892

10

Columns 7 through 12

0.0164

0.0125

0.0099

0.0080

0.0066

0.0056

1.5785

1.9731

2.3677

3.9462

4.3408

4.7354

Column 13

0.0047

R2 =

1.0e+08 *

Columns 1 through 6

0.3946

0.7892

1.1839

Columns 7 through 12

2.7623

3.1569

Column 13

5.1300

Ratio3 =

1.0e-21 *

3.5516

11

Columns 1 through 6

0.2034

0.0254

0.0075

0.0032

0.0016

0.0009

0.0002

0.0002

0.0001

Columns 7 through 12

0.0006

0.0004

0.0003

Column 13

0.0001

Codes in Command Window: plot(velocity,Ratio1); hold on plot(velocity,Ratio2); hold on plot(velocity,Ratio3)

Graph Obtained:

12

3. By volume, dry air contains 78.09% nitrogen, 20.95% oxygen, 0.93% argon, 0.04% carbon dioxide, and small amounts of other gases. Express this Data in terms of a Pie Chart for better understanding. MATLAB Codes: X = [78.09 0.93 20.95 0.04]; >> labels = {'Nitrogen' 'Argon' 'Oxygen' 'Carbondioxide and Other Gases'}; >> pie(X,labels)

Result:

Interpretation: From the above graphs, we not only get an idea of variation of different parameters with respect to a constantly varying variable but also get to compute the value of required parameters at any given point without mathematical interpolation.

13

4. Various Flight paths have been given along with their magnitudes of ranges it travelled. Interpret the data using visualisation.

MATLAB Code: z = eig(randn(10)); compass(z)

Result:

Interpretation: Thus the Flight Paths of Various Aircrafts taking off an airport have been visualized and interpreted.

14

5. The error occurred in measuring the airspeed by the onboard Pitot-Static tube is calculated for different speeds in terms of Mach Number. Interpret the data in MATLAB.

MATLAB Code: x = 0:0.1:1; errorbar(x,exp(–x), 0.5*rand(1,length(x)),’d’)

Result:

Interpretation: Hence the error obtained in calculation of Airspeed for different Mach Numbers have been interpreted and analysed.

15

3D Plot Theory: Some types of 3D graphical interpretation and its examples are as follows:

Mesh Plot: The mesh function creates a wireframe mesh. By default, the color of the mesh is proportional to the surface height. z = peaks(25); figure mesh(z) Surface Plot: The surf function is used to create a 3-D surface plot. surf(z) colormap(jet)

Surface Plot (with Shading): The surfl function creates a surface plot with colormap-based lighting. For smoother color transitions, use a colormap with linear intensity variation such as pink. surfl(z) colormap(pink) % change color map shading interp

% interpolate colors across lines and faces

Contour Plot: The contour function is used to create a plot with contour lines of constant value. contour(z,16) colormap default

% change color map

Quiver Plot: The quiver function plots 2-D vectors as arrows. x = -2:.2:2; y = -1:.2:1;

16

[xx,yy] = meshgrid(x,y); zz = xx.*exp(-xx.^2-yy.^2); [px,py] = gradient(zz,.2,.2);

quiver(x,y,px,py) xlim([-2.5 2.5])

Problems:

1. The Variation of Heat of the jet flow coming out of the Nozzle of a turbojet engine is given by the equation: 5 𝑄= 1 + 𝑇 2 + 𝑇𝑖2 Where Q is the heat generated, T is the nozzle exit temperature and Ti is the Nozzle inlet temperature. Plot the variation of Heat with these two parameters.

MATLAB CODES: X = -3 : .1 : 3; [x,y] = meshgrid(X,X); z = 1./(3+x.^2+y.^2); surfl(z) shading interp colormap hot xlabel('x'); ylabel('y');

17

GRAPH:

Interpretation: From this Graph of variation of Heat generated to Temperatures at inlet and outlet, we get to know the locations of Hotspots throughout the flow in nozzle at high speeds of the turbojet aircrafts.

2. Visualise the Earth’s Surface using the 3D graph along with the Latitude and Longitude of Earth in the graph.

MATLAB CODE: sphere(50) axis 'equal'

18

GRAPH:

Interpretation: From this visualisation, we get a rough 3D image of Earth which can be used to plot various points such as the poles, magnetic poles, etc. Also they can be used to track satellite trajectories at the time of launch etc.

3. The variation of Dynamic Viscosity of Ethanol flowing at 25ᵒ inside a tube is given by the relation: µ = cos 4𝜋𝜌 + 2 where ρ is the density of the ethanol in the flow. Plot the variation of density and Dynamic Viscosity and analyse.

MATLAB Codes: z = 0: .03 : 1; r = cos(4*pi*z)+2; cylinder(r)

19

GRAPH:

Interpretation: As the dynamic viscosity varies with the density with a function of Cosine. We obtain the graph with varies the same way harmonically again and again throughout the flow. It is a point to be notes that density does not remains constant and depends on several factors such as the temperature, pressure etc.

20

4. When air passes through a series of Expansion waves, there is change in temperature along with various parameters related to them. Take random values of Two expansion waves and use a suitable plot to depict the relation between them.

MATLAB CODE: x = [0 2.5; 5 2.5; 5 2.5; 0 2.5]; y = [0 0; 0 -1; 0 -1; 0 0]; z = [0 0; 0 0; 2 2; 2 2]; fill3(x,y,z, rand(4,2)) xlabel('x'); ylabel('y'); zlabel('z'); view(120, 50) grid

GRAPH:

21

5. The Infrared sensors of combat Heli-carriers used advanced algorithms to detect an obstacle such as an elevated mountain, transforms it into an polynomial equations and interprets the height and other parameters for safer cruise. Consider a equation of the mountain depicted by the algorithm to be: 1 𝑧= log 8 + 𝑥 3 + 3𝑦 2 Plot the 3D visualisation of this obstacle as done by the On-board Systems.

MATLAB CODE: X = -3:.1:3; [x,y] = meshgrid(X,X); z = 1./(8+x.^2+3*y.^2); contour3(z) xlabel('x'); ylabel('y');

GRAPH:

Interpretation: The Heli-carrier’s on board computers uses this data to find the altitude of the obstacle and find an alternate safe Cruise path for it to fly.

22

Eigen Values and Vectors

Theory: e = eig(A) returns a column vector containing the eigenvalues of square matrix A. [V,D] = eig(A) returns diagonal matrix D of eigenvalues and matrix V whose columns are the corresponding right eigenvectors, so that A*V = V*D. [V,D,W] = eig(A) also returns full matrix W whose columns are the corresponding left eigenvectors, so that W'*A = D*W'. The eigenvalue problem is to determine the solution to the equation Av = λv, where A is an nby-n matrix, v is a column vector of length n, and λ is a scalar. The values of λ that satisfy the equation are the eigenvalues. The corresponding values of v that satisfy the equation are the right eigenvectors. The left eigenvectors, w, satisfy the equation w’A = λw’. e = eig(A,B) returns a column vector containing the generalized eigenvalues of square matrices A and B. [V,D] = eig(A,B) returns diagonal matrix D of generalized eigenvalues and full matrix V whose columns are the corresponding right eigenvectors, so that A*V = B*V*D. [V,D,W] = eig(A,B) also returns full matrix W whose columns are the corresponding left eigenvectors, so that W'*A = D*W'*B. The generalized eigenvalue problem is to determine the solution to the equation Av = λBv, where A and B are n-by-n matrices, v is a column vector of length n, and λ is a scalar. The values of λ that satisfy the equation are the generalized eigenvalues. The corresponding values of v are the generalized right eigenvectors. The left eigenvectors, w, satisfy the equation w’A = λw’B. [___] = eig(A,balanceOption), where balanceOption is 'nobalance', disables the preliminary balancing step in the algorithm. The default for balanceOption is 'balance', which enables balancing. The eig function can return any of the output arguments in previous syntaxes. [___] = eig(A,B,algorithm), where algorithm is 'chol', uses the Cholesky factorization of B to compute the generalized eigenvalues. The default for algorithm depends on the properties of A and B, but is generally 'qz', which uses the QZ algorithm. If A is Hermitian and B is Hermitian positive definite, then the default for algorithm is 'chol'. [___] = eig(___,eigvalOption) returns the eigenvalues in the form specified by eigvalOption using any of the input or output arguments in previous syntaxes. Specify eigvalOption as 'vector' to return the eigenvalues in a column vector or as 'matrix' to return the eigenvalues in a diagonal matrix.

23

Problems: 1. The Stress on X,Y and Z axes of a given material Niantic-75 used in the manufacturing of Combustion Chamber is given by the Eigen Values of this Matrix. Compute using MATLAB.

MATLAB Code: MatrixA=[3 2 4;2 0 2;4 2 3] e=eig(MatrixA)

Result: MatrixA =

3

2

4

2

0

2

4

2

3

e= -1.0000 -0.4721 8.4721

Interpretation: The Stress Values are 1, 0.4721 and 8.4721 N/m2 on the 3-Axis respectively. The Negative Sign denotes the Nature of the stress which is either Compressive or Tensile.

24

2. A single Crystal of Diamond used for cutting Aircraft Parts is loaded as follows on its [100] , [010] and [001] axes. What are the normal stresses on the orthogonal [111], [110], [112] axes?

MATLAB Code: MatrixB=[100 20 0;20 0 20;0 20 100]; >> e=eig(MatrixB)

Result: e= -7.4456 100.0000 107.4456

Interpretation: The Normal Stress along the orthogonal plane of [111], [110], [112] are -7.4456, 100.0, 107.4456 respectively.

25

Polynomial Evaluation: Theory: y = polyval(p,x) evaluates the polynomial p at each point in x. The argument p is a vector of length n+1 whose elements are the coefficients (in descending powers) of an nth-degree

polynomial: p(x)=p1xn+p2xn−1+...+pnx+pn+1. The polynomial coefficients in p can be calculated for different purposes by functions like polyint, polyder, and polyfit, but you can specify any vector for the coefficients. To evaluate a polynomial in a matrix sense, use polyvalm instead. [y,delta] = polyval(p,x,S) uses the optional output structure S produced by polyfit to generate error estimates. delta is an estimate of the standard error in predicting a future observation at x by p(x).

y = polyval(p,x,[],mu) or [y,delta] = polyval(p,x,S,mu) use the optional output mu produced by polyfit to center and scale the data. mu(1) is mean(x), and mu(2) is std(x). Using these values, polyvalcenters x at zero and scales it to have unit

standard deviation. This centering and scaling transformation improves the numerical properties of the polynomial. We use Polyint function in orde to perform the Integration Functions in MATLAB.

Problems: 1. The weight of an ideal round cut diamond which is used in cutting the Aluminium for Aircraft can be modelled by the equation: w=0.0021d3-0.060d2+0.23d where w is the weight of the diamond in carats and d is the diameter in millimetres. For Cutting a specific Part, diamond of 30 Millimetres is required. Determine the Weight of diamond required.

26

MATLAB Codes: p=[0.0021 -0.060 0.023]; x=[30]; y=polyval(p,x)

Result: y =0.1130

Interpretation: This Result gives is the Value of the Volume of the Diamond at x=30 mm. Hence a diamond with volume of 0.1130 mm3 is required for the specific cutting process.

2. Consider a Symmetric Airfoil whose Mid-section is in the form of the rectangle with sides of magnitude of Z, (Z+8), (Z+2) respectively. The volume is given by the equation: V=4z3+16z2+12z+15 A particular NACA Series Airfoil should have X=12 mm of dimension. Find the Volume of the Rectangle used in the Airfoil.

MATLAB Codes: Volume=[4 16 12 15]; z=[12]; y=polyval(Volume,z)

Result: y = 9375

Interpretation: Hence, The Volume of the rectangular section of the required airfoil is 9375 mm3.

27

CURVE FITTING: Theory: To interactively fit a curve, follow the steps in this simple example: 1. Load some data at the MATLAB® command line. load hahn1 2. Open the Curve Fitting app. Enter: cftool 3. In the Curve Fitting app, select X Data and Y Data. Curve Fitting app creates a default interpolation fit to the data. 4. Choose a different model type using the fit category drop-down list, e.g., select Polynomial. 5. Try different fit options for your chosen model type. 6. Select File > Generate Code. Curve Fitting app creates a file in the Editor containing MATLAB code to recreate all fits and plots in your interactive session. Alternatively, Polynomial curve fitting Syntax p = polyfit(x,y,n) [p,S] = polyfit(x,y,n) [p,S,mu] = polyfit(x,y,n) Description: p = polyfit(x,y,n)

returns the coefficients for a polynomial p(x) of degree n that is a best fit (in a least-squares sense) for the data in y. The coefficients in p are in descending powers, and the length of p is n+1 p(x)=p1xn+p2xn−1+...+pnx+pn+1. [p,S] = polyfit(x,y,n) also returns a structure S that can be used as an input to polyval to obtain error estimates. [p,S,mu] = polyfit(x,y,n) also returns mu, which is a two-element vector with and scaling values. mu(1) is mean(x), and mu(2) is std(x). Using these values, polyfit centers x at zero and scales it to have unit standard deviation,

ˆx= ‾x x . x− σ

centering

28

This centring and scaling transformation improves the numerical properties of both the polynomial and the fitting algorithm.

Questions: 1. The Coordinates of a part of the Horizontal Stabilizer of an Aircraft is given as follows: X Y 1

14

2

27

3

40

4

55

5

68

6

75

7

84

8

91

9

104

10

125

11

136

Plot the Graph using Curve Fitting Tool.

MATLAB Codes: X_coordinates = [1 2 3 4 5 6 7 8 9 10 11]; >> Y_coordinates = [14 27 40 55 68 75 84 91 104 125 136]; >> polyfit(X_coordinates,Y_coordinates,1)

29

Result: ans =

11.6545

4.5273

Interpretation: From the graph, We come to know that the plot of the given coordinates is of the Horizontal Stabilizer of the given aircraft.

2. The Variation between the Degree Shift of the Flaps mounted on the Wings of the aircraft which is used to increase the camber of the Aircraft wing with the Percentage increase in Coefficient of Lift varies with each other in a Parabolic form. With the given Data, prove it graphically. X Y 0

1.8

1

1.98

2

1.9926

3

2.07

4

2.67

30

5

2.88

6

3.00

MATLAB Codes: X_coordinates = [0 1 2 3 4 5 6]; >> Y_coordinates = [1.8 1.98 1.9926 2.07 2.67 2.88 3.00]; polyfit(X_coordinates,Y_coordinates,1)

Result: ans = 0.0206

0.0933

1.7938

31

Solving Algebraic and Transcendental Equations: Theory: If eqn is an equation, solve(eqn, x) solves eqn for the symbolic variable x. Use the == operator to specify the familiar quadratic equation and solve it using solve. syms a b c x eqn = a*x^2 + b*x + c == 0; solx = solve(eqn, x) solx = -(b + (b^2 - 4*a*c)^(1/2))/(2*a) -(b - (b^2 - 4*a*c)^(1/2))/(2*a) solx is a symbolic vector containing the two solutions of the quadratic equation. If the input eqn is an expression and not an equation, solve solves the equation eqn == 0.

Questions: 1. The Error in Calculating the Pressure at a Given Altitude is given by the equation: 8x3+12x2+5x+2=0 Find the value of X. As 3 different values are obtained, the highest possible error Magnitude is considered.

MATLAB Code: Polynomial=[8 12 5 2]; >> root=roots(Polynomial)

Result: root = -1.1448 + 0.0000i -0.1776 + 0.4322i -0.1776 - 0.4322i

Interpretation: Hence, the three required roots of the given polynomial are obtained.

32

2. The Given quadratic equation is the equation obtained while solving for a constant that determines the Oswald’s Efficiency value which is used in the calculation of various parameters in Flight Dynamics. Solve it using MATLAB. 12x2+4x+16=0

MATLAB Code: Polynomial=[12 4 16]; >> root=roots(Polynomial)

Result: root = -0.1667 + 1.1426i -0.1667 - 1.1426i

Interpretation: We get the 2 different values of Oswald’s Constant from the MATLAB.

Questions on Transcendental Equations: 1. The Resonance occurring on a Aircraft Wing on Flight is given by the equation: 𝑓(𝛽) = cosh 𝛽 ∗ cos 𝛽 + 1 = 0 Where 𝛽𝑖4 = (2𝜋𝑓𝑖 )2 (

𝑚𝐿3 ) 𝐸𝐼

Where 𝑓𝑖 = ith natural frequency m= mass of the beam L= Length of the beam E= Modulus of Elasticity I=Moment of inertia of the Cross-section Find the Angles at which resonance occurs at the wing of this Airplane.

33

MATLAB Code: >> syms x >> Eqn2 = cosh(x)*cos(x)+1 == 0; >> solx = solve(Eqn2,x)

Result: solx = -212.05750411731104359622842837137

Interpretation: The given value is close to 1.177 Radians at which there can be observed radians on the wing of this Aircraft.

2. The Specific Fuel input to the Engine for an Commercial Aircraft is given by the Lambert’s Equation: 𝑎𝑥 + 𝑏 log 𝑥 + 𝑐 √𝑥 + 𝑑 = 0 For a particular Aircraft Flying at 30,000 ft, the Specific Fuel Input is given by the equation: 3𝑥 + 5 log 𝑥 + 9√𝑥 + 10 = 0 Solve the equation using MATLAB.

MATLAB Code: Eqn = 3*y+5*log(y)+9*sqrt(y)+10 ==0; >> solx=solve(Eqn,y)

Result: solx =0.078092206567782850204281665495309

Interpretation: Hence, The Specific Fuel Input for this Aircraft at 30,000 ft is 0.078092.

34

System of Linear Equations Theory:

Questions: 1. The Three-Moment Equation for determining the Stress across a Beam is given by: 6𝑎1 𝑥̅1 6𝑎2 𝑥̅2 𝑀𝑎 𝐿1 + 2𝑀𝐵 (𝐿1 + 𝐿2 ) + 𝑀𝐶 𝐿2 = + 𝐿1 𝐿2

The equations obtained after substituting the values are: 1. 22MA + 5MB=172.36 2. 5MA + 18MB=95.32 Solve for both the Moments.

35

MATLAB CODE: MAtA=[22 5;5 18]; MatB=[172036;95.32] X=linsolve(MatA,MatB)

Solution: X= 7.0778 3.3295

Interpretation: Hence we obtain the Values of MA and MB to be 7.0778 and 3.3295 respectively.

2. The Three-Moment Equation for determining the Stress across a Beam is given by: 6𝑎1 𝑥̅1 6𝑎2 𝑥̅2 𝑀𝑎 𝐿1 + 2𝑀𝐵 (𝐿1 + 𝐿2 ) + 𝑀𝐶 𝐿2 = + 𝐿1 𝐿2 The equations obtained after substituting the values are: 1. 24MA + 16MB + 3MC= 62.5 2. 5MA + 7MB + 12MC= 158.63 3. 8MA + 2MB + 9Mc= 52.32

MATLAB Codes: MAtA=[24 16 3;5 7 12;8 2 9]; MatB=[62.5;158.83;52.32] X=linsolve(MatA,MatB)

36

Result: X= -6.5268 12.0195 8.9440

Interpretation: Hence we obtain the Values of MA, MB and MC to be 6.5268, 12.0944 and 8.944 respectively.

3. The mass flow rate of air be 2 kg/s and 5 kg/s and mass flow rate of fuel be 0.001kg/s and 0.002kg/s and the thrust produced be 12,487N and 14,698N. Calculate the mass flow rate of air and fuel for maximum thrust generated by the Aircraft.

MATLAB Code: a=[ 2 5;0.001 0.002]; b=[12,487; 14,698]; linsolve(a,b)

Results: ans = 3.56212 0.01245

Interpretation: Hence the Mass flow rate and fuel required for generating maximum thrust has been calculated.

37

Numerical Integration of Single/Double/Triple Integral Theory: Single Integration: Syntax q = integral(fun,xmin,xmax) q = integral(fun,xmin,xmax,Name,Value)

Description q = integral(fun,xmin,xmax) numerically integrates function fun from xmin to xmax using global adaptive quadrature and default error tolerances. q = integral(fun,xmin,xmax,Name,Value) specifies additional options with one or more Name,Value pair arguments. For example, specify 'WayPoints' followed by a vector of real or complex numbers to indicate specific points for the integrator to use.

Double Integration: Syntax q = integral2(fun,xmin,xmax,ymin,ymax) q = integral2(fun,xmin,xmax,ymin,ymax,Name,Value)

Description q = integral2(fun,xmin,xmax,ymin,ymax) approximates the integral of the function z = fun(x,y) over the planar region xmin ≤ x ≤ xmax and ymin(x) ≤ y ≤ ymax(x). q = integral2(fun,xmin,xmax,ymin,ymax,Name,Value) specifies additional options with one or more Name,Value pair arguments.

Triple Integration: Syntax q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax) q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax,Name,Value)

38

Description q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax) approximates the integral of the function z = fun(x,y,z) over the region xmin ≤ x ≤ xmax, ymin(x) ≤ y ≤ ymax(x) and zmin(x,y) ≤ z ≤ zmax(x,y). q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax,Name,Value) specifies additional options with one or more Name,Value pair arguments.

Questions: 1. The Spar Extension of the Airfoil NACA 2414 in a section middle of the Airfoil is evaluated by integrating the Equation of the Curvature by keeping the limits as distances of both the ends from the Root Chord Locations respectively. Get the Spar Extension value of a section whose endpoints are at a distances of 4 and 9 meters respectively from the Root Chord respectively.

MATLAB CODES: >> fun = @(x) (1+x)./(x+2.*sqrt(x)-3); >> q=integral(fun,4,9)

Result: q=

4.4280

Interpretation: From the Integration, we get to know that a spar of length 4.428m is required to be placed at that area of the airfoil.

39

2. A satellite is launched and is in a circular orbit at an altitude of 36000 kms in a geostationary orbit. The satellite has to transferred to a elliptical orbit using Hoffmann’s Transfer process. For this, the onboard impulse Thrusters2*cos(x) have to generate impulse thrust to rotate the satellite at a particular angle. In Hoffmann’s transfer, the satellite has to rotated to an angle of 90 degrees. The amount of thrust required is obtained by integrating the equation with the desired limits.

MATLAB Codes: >> fun = @(x) sin(x).*4.*cos(x).*exp(2*cos(x)+1); >> r=integral(fun,0,pi/2)

Result: r= 22.8038

Interpretation: Hence an Impulse force of 22.80 Newtons has to be propelled by the onboard thrusters in order to switch to the desired orbit. 3. When an Aircraft flies in supersonic speeds, there is heat generation on the outer surfaces of the aircraft. Considering the surface to be a 2 dimensional object, The heat generation is taken in 2 components which is in X and Y direction. On a particular areal section of an aircraft flying at Mach 2.8, The heat generation in Kilojoules is given by integrating the given equation with the area magnitudes as limits. Obtain the heat generated on this Particular aircraft.

MATLAB Code: >> func = @(x,y) 1+8.*x.*y; >> s=integral2(func,0,3,1,2)

Result: s = 57

40

Interpretation: The Heat generated on that particular part of the surface is obtained to be 57 Kilojoules by the method of double Integration. 4. Consider a SR-71, a supersonic Aircraft flying at an altitude of 50000 ft above sea level at a speed of Mach 3.9. A shock wave is formed in the form of a 3 dimensional cone. The volume inside that Shock Cone is known as the zone of Action where there is disturbance of the air molecules in it. Integrating the below equation gives the Volume of the Mach Cone which is termed as the zone of action in this case of supersonic flight. Integrate and obtain the Volume of zone of action for this case.

MATLAB Code: >> funct = @(x,y,z) x.*y.*z; >> xmin=0 >> xmax=1; >> ymin=0; >> ymax=1; >> zmin=@(x,y) sqrt(x.^2+y.^2); >> zmax= 2; >> q = integral3(funct,xmin,xmax,ymin,ymax,zmin,zmax)

Result: q= 0.375

Interpretation: The Volume of the Mach cone or the zone of action is determined to be 0.375 Cubic kilometres.

41

Integration of Numeric Data Theory: Trapz: Syntax Q = trapz(Y) Q = trapz(X,Y) Q = trapz(___,dim)

Description Q = trapz(Y) computes the approximate integral of Y via the trapezoidal method with unit spacing. The size of Y determines the dimension to integrate along: • • •

• •

If Y is a vector, then trapz(Y) is the approximate integral of Y. If Y is a matrix, then trapz(Y) integrates over each column and returns a row vector of integration values. If Y is a multidimensional array, then trapz(Y) integrates over the first dimension whose size does not equal 1. The size of this dimension becomes 1, and the sizes of other dimensions remain unchanged. Q = trapz(X,Y) integrates Y with respect to the coordinates or scalar spacing specified by X. If X is a vector of coordinates, then length(X) must be equal to the size of the first dimension of Y whose size does not equal 1. If X is a scalar spacing, then trapz(X,Y) is equivalent to X*trapz(Y). Q = trapz(___,dim) integrates along the dimension dim using any of the previous syntaxes. You must specify Y, and optionally can specify X. If you specify X, then it can be a scalar or a vector with length equal to size(Y,dim). For example, if Y is a matrix, then trapz(X,Y,2) integrates each row of Y.

Cumtrapz: Syntax Q = cumtrapz(Y) Q = cumtrapz(X,Y)

42

Q = cumtrapz(___,dim)

Description

• • •

• •

Q = cumtrapz(Y) computes the approximate cumulative integral of Y via the trapezoidal method with unit spacing. The size of Y determines the dimension to integrate along: If Y is a vector, then cumtrapz(Y) is the cumulative integral of Y. If Y is a matrix, then cumtrapz(Y) is the cumulative integral over each column. If Y is a multidimensional array, then cumtrapz(Y) integrates over the first dimension whose size does not equal 1. Q = cumtrapz(X,Y) integrates Y with respect to the coordinates or scalar spacing specified by X. If X is a vector of coordinates, then length(X) must be equal to the size of the first dimension of Y whose size does not equal 1. If X is a scalar spacing, then cumtrapz(X,Y) is equivalent to X*cumtrapz(Y). Q = cumtrapz(___,dim) integrates along the dimension dim using any of the previous syntaxes. You must specify Y, and optionally can specify X. If you specify X, then it can be a scalar or a vector with length equal to size(Y,dim). For example, if Y is a matrix, then cumtrapz(X,Y,2) cumulatively integrates each row of Y.

SQuestions: 1. In the experimentation process of measure surface pressure distribution across a wing surface placed in a supersonic wind tunnel, the Net Lift generated is obtained by integrating the a depending pressure factor with respect to the distance or the position of the sensor with respect to the leading edge tip with the help on integration with unit spacing. The data of the dependent factor and the sensor position is given as follows: Location of probes 0 0.1 0.2 0.4 0.8 1.5 2.7 3.9 5.2 6.5 7.8 9

Depending Pressure Factor 1 0.6111 0.3889 0.222 0.556 0.333 0.7222 0.555 0.3889 0.6578 0.5642 0.1203

43

10

0.892

Find the value of lift generated using the given data with the help of Unit Spacing Integration.

MATLAB Code: >> location = [0 0.1 0.2 0.4 0.8 1.5 2.7 3.9 5.2 6.5 7.8 9 10] >> Factor = [1 0.6111 0.3889 0.222 0.556 0.333 0.7222 0.555 0.3889 0.6578 0.5642 0.1203 0.892] >> y=trapz(location,Factor)

Result: y=

5.0629

Interpretation: Hence a net Lift of 5 Newtons is generated on the scaled-down wing model which is being tested in the supersonic Wind Tunnel.

2. The braking system of the landing Wheels of a Long haul Boeing 787 is being tested. The Static Braking Distance of the wheels is dependent on the Coefficient of Static Friction between the disk brakes used in the flight. The Total braking distance is obtained by integrating the Coefficient of Static Friction values per unit time in minutes. Gives the values of Coefficient of Static friction of a trial. Compute the Static Braking Distance of the flight.

Time(in minutes)

Coeff. Of Static Friction

1

1.13

2

1.68

3

1.86

4

2.032

5

2.521

6

3.12

44

7

4.213

8

5.658

9

7.1453

10

9.546

MATLAB Code: >> Friction_factor = [1.13 1.68 1.86 2.032 2.521 3.12 4.213 5.658 7.1453 9.546]; >> x=cumtrapz(Friction_factor)

Result: x=

Columns 1 through 6

0

1.4050

3.1750

5.1210

7.3975 10.2180

Columns 7 through 10

13.8845 18.8200 25.2217 33.5673

Interpretation: The Static Braking Distances for each time Interval in minutes have been computed for this Flight.

45

Solving Ordinary (Linear) Differential Equations Theory: [t,y] = ode45(odefun,tspan,y0), where tspan = [t0 tf], integrates the system of differential equations y'=f(t,y) from t0 to tf with initial conditions y0. Each row in the solution array y corresponds to a value returned in column vector t. All MATLAB® ODE solvers can solve systems of equations of the form y'=f(t,y), or problems that involve a mass matrix, M(t,y)y'=f(t,y). The solvers all use similar syntaxes. The ode23s solver only can solve problems with a mass matrix if the mass matrix is constant. ode15s and ode23t can solve problems with a mass matrix that is singular, known as differential-algebraic equations (DAEs). Specify the mass matrix using the Mass option of odeset. ode45 is a versatile ODE solver and is the first solver you should try for most problems. However, if the problem is stiff or requires high accuracy, then there are other ODE solvers that might be better suited to the problem. See Choose an ODE Solver for more information. [t,y] = ode45(odefun,tspan,y0,options) also uses the integration settings defined by options, which is an argument created using the odeset function. For example, use the AbsTol and RelTol options to specify absolute and relative error tolerances, or the Mass option to provide a mass matrix. [t,y,te,ye,ie] = ode45(odefun,tspan,y0,options) additionally finds where functions of (t,y), called event functions, are zero. In the output, te is the time of the event, ye is the solution at the time of the event, and ie is the index of the triggered event. For each event function, specify whether the integration is to terminate at a zero and whether the direction of the zero crossing matters. Do this by setting the 'Events'property to a function, such as myEventFcn or @myEventFcn, and creating a corresponding function: [value,isterminal,direction] = myEventFcn(t,y). For more information, see ODE Event Location. sol = ode45(___) returns a structure that you can use with deval to evaluate the solution at any point on the interval [t0 tf]. You can use any of the input argument combinations in previous syntaxes. [t,y] = ode23(odefun,tspan,y0), where tspan = [t0 tf], integrates the system of differential equations y'=f(t,y) from t0 to tf with initial conditions y0. Each row in the solution array y corresponds to a value returned in column vector t. All MATLAB® ODE solvers can solve systems of equations of the form y'=f(t,y), or problems that involve a mass matrix, M(t,y)y'=f(t,y). The solvers all use similar syntaxes. The ode23s solver only can solve problems with a mass matrix if the mass matrix is constant. ode15s and ode23t can solve problems with a mass matrix that is singular, known as differential-algebraic equations (DAEs). Specify the mass matrix using the Mass option of odeset.

46

[t,y] = ode23(odefun,tspan,y0,options) also uses the integration settings defined by options, which is an argument created using the odeset function. For example, use the AbsTol and RelTol options to specify absolute and relative error tolerances, or the Mass option to provide a mass matrix. [t,y,te,ye,ie] = ode23(odefun,tspan,y0,options) additionally finds where functions of (t,y), called event functions, are zero. In the output, te is the time of the event, ye is the solution at the time of the event, and ie is the index of the triggered event. For each event function, specify whether the integration is to terminate at a zero and whether the direction of the zero crossing matters. Do this by setting the 'Events'property to a function, such as myEventFcn or @myEventFcn, and creating a corresponding function: [value,isterminal,direction] = myEventFcn(t,y). For more information, see ODE Event Location. sol = ode23(___) returns a structure that you can use with deval to evaluate the solution at any point on the interval [t0 tf]. You can use any of the input argument combinations in previous syntaxes.

Questions: 1. Calculate the rate of temperature change with respect to time in the conduction heat transfer of rod of length 10 metre. The material has thermal conductivity of 100W/m2 K and the area of the rod be 12*10-2 metre2.

MATLAB Code: Q=(10:1:100) k=100 A=12*10^-2 t=(Q./-k*A) tspan=[0 400] y0=0 [Q,t] = ode45(@(Q,t)Q./-k*A, tspan, y0) [Q,t]=meshgrid(-10:1:10) R=Q./t surf(Q,t,R) xlabel('Heat flux') ylabel('Temperature change') zlabel('resistance') title('Variation of heat transfer by conduction wrt heat flux, temperature,resistance along the rod')

47

Results:

2. There is spherical tank for storing water. The tank is filled through a hole in the top and drained through a hole in the bottom. If the tank’s radius is r, you can use integration to show that the volume of water in the tank as a function of its height h is given by 𝑉(ℎ) = 𝜋ℎ2 𝑟 − 𝜋

ℎ3 3

Find the Draining rate from the tank.

Solution: Torricelli’s principle states that the liquid ow rate through the hole is proportional to the square root of the height h. Further studies in fluid mechanics have identified the relation more precisely, and the result is that the volume ow rate through the hole is given by 𝑞 = 𝐶𝑑 𝐴√2𝑔ℎ where A is the area of the hole, g is the acceleration due to gravity, and Cd is an experimentally determined value that depends partly on the type of liquid. For water, Cd=0.6 is a common value. We can use the principle of conservation of mass to obtain a

48

differential equation for the height h. Applied to this tank, the principle says that the rate of change of liquid volume in the tank must equal the ow rate out of the tank; that is,

Draining of a spherical tank

𝑑ℎ 0.0334√ℎ =− 𝑑𝑡 10ℎ − ℎ2

MATLAB Code: t=(0:1:247); h=(0:1:10); [t, h]= ode45 (@(t,h) (sqrt(h))/(10*h+h^2),[0,2475], 10) plot(t,h) ,xlabel('Time (sec)'), ylabel('Height (ft)') [t, y] = ode45(@(t,y) 2*exp(-10*t) ,tspan, 2);

Results:

49

Interpretation: The resulting plot is shown in Figure. Note how the height changes more rapidly when the tank is nearly full or nearly empty. This is to be expected because of the effects of the tank’s curvature. The tank empties in 2475 sec, or 41 min. This value is not grossly different from our rough estimate of 52 min, so we should feel comfortable accepting the numerical results. The value of the nominal time of 2475 sec was found by increasing the nominal time until the plot showed that the height became 0.

3. The pendulum shown in Figure 9.4–1 consists of a concentrated mass m attached to a rod whose mass is small compared to m. The rod’s length is L. The equation of motion for this pendulum is which is linear and has the solution

𝑔 𝜃(𝑡) = 𝜃(0) cos √ 𝑡 𝐿

MATLAB Code: function xdot = pendulum(t,x) g = 9.81; L = 1; xdot = [x(2); -(g/L)*sin(x(1))]; [ta, xa] = ode45(@pendulum, [0,5], [0.5, 0]); [tb, xb] = ode45(@pendulum, [0,5], [0.8*pi, 0]); plot(ta, xa(:,1), tb,xb(:,1)), xlabel ('Time (s)') ylabel('Angle (rad)'), gtext('Case 1'), gtext('Case 2')

Results:

50

Solving Ordinary (Nonlinear) Differential Equations Theory: In mathematics and science, a nonlinear system is a system in which the change of the output is not proportional to the change of the input. Nonlinear problems are of interest to engineers, biologists, physicists, mathematicians, and many other scientists because most systems are inherently nonlinear in nature. Nonlinear dynamical systems, describing changes in variables over time, may appear chaotic, unpredictable, or counterintuitive, contrasting with much simpler linear systems. Typically, the behaviour of a nonlinear system is described in mathematics by a nonlinear system of equations, which is a set of simultaneous equations in which the unknowns (or the unknown functions in the case of differential equations) appear as variables of a polynomial of degree higher than one or in the argument of a function which is not a polynomial of degree one. In other words, in a nonlinear system of equations, the equation(s) to be solved cannot be written as a linear combination of the unknown variables or functions that appear in them. Systems can be defined as nonlinear, regardless of whether known linear functions appear in the equations. In particular, a differential equation is linear if it is linear in terms of the unknown function and its derivatives, even if nonlinear in terms of the other variables appearing in it. As nonlinear dynamical equations are difficult to solve, nonlinear systems are commonly approximated by linear equations (linearization). This works well up to some accuracy and some range for the input values, but some interesting phenomena such as solitons, chaos, and singularities are hidden by linearization. It follows that some aspects of the dynamic behaviour of a nonlinear system can appear to be counterintuitive, unpredictable or even chaotic. Although such chaotic behaviour may resemble random behaviour, it is in fact not random. For example, some aspects of the weather are seen to be chaotic, where simple changes in one part of the system produce complex effects throughout. This nonlinearity is one of the reasons why accurate long-term forecasts are impossible with current technology. NONSTIFF ODES This page contains two examples of solving non-stiff ordinary differential equations using ode45. MATLAB has three solvers for non-stiff ODEs. •

ode45

•

ode23

•

ode113

For most non-stiff problems, ode45 performs best. However, ode23 is recommended for problems that permit a slightly cruder error tolerance or in the presence of moderate stiffness. Likewise, ode113 can be more efficient than ode45 for problems with stringent error tolerances.

51

If the non-stiff solvers take a long time to solve the problem or consistently fail the integration, then the problem might be stiff. See for more information. STIFF ODES This page contains two examples of solving stiff ordinary differential equations using ode15s. MATLAB® has four solvers designed for stiff ODEs. •

ode15s

•

ode23s

•

ode23t

•

ode23tb

For most stiff problems, ode15s performs best. However, ode23s, ode23t, and ode23tb can be more efficient if the problem permits a crude error tolerance.

Questions: 1. The rate of fuel consumption of a turbojet engine varies constantly with respect to amount of thrust generated by the Engine. The Fuel inlet valve has to constantly controlled to control the fuel input. The Fuel consumption varies with the thrust differentially with the equation: 𝑦 = 𝑥′2 − 3𝑥 + 27

Solve the equation for data of Thrust and Fuel input Values.

52

MATLAB Code: Fun=@(x,y) y=x1^2-3*x+27 A= [0 12]; B = 27; [x,y] = ode45(@(x,y) (fun,A, B) plot(x,y,'-o')

Result: A=

0.121 1.3721 1.7443 1.1164 1.4885 3.3492 5.2098 7.0705 8.9312 12.1312 15.3312

53

18.5312 21.7312 24.9312 28.1312 31.3312 34.5312 37.7312 40.9312 44.1312 47.3312 50.5312 53.7312 56.9312 60.1312 63.3312 66.5312 69.7312 72.9312 76.1312 79.3312 82.5312 85.7312

54

88.9312 92.1312 95.3312 98.5312 101.7312 104.9312 108.1312 111.3312 114.5312 117.7312 120.9312 124.1312 125.0984 126.0656 127.0328 128.0000

B=

1.0+06 *

55

0.0200 0.0211 0.0223 0.0237 0.0252 0.0353 0.0491 0.0668 0.0884 0.1344 0.1918 0.2606 0.3408 0.4323 0.5352 0.6494 0.7750 0.9120 1.0603 1.2200 1.3911 1.5736

56

1.7674 1.9726 2.1891 2.4170 2.6563 2.9069 3.1689 3.4423 3.7270 4.0232 4.3306 4.6495 4.9797 5.3212 5.6742 6.0385 6.4142 6.8012 7.1996 7.6094 8.0305 8.4630

57

8.9069 9.0433 9.1807 9.3192 9.4587

Graph:

Interpretation: Thus different values of Fuel input factor for different thrusts generated by the engines are being computed and plotted for observations.

58

2. The Icing is a common phenomenon occurring on aircraft body flying at high altitudes of at cold temperatures. Ina scenario, where Gulfstream 4 is flying at higher altitudes, the ice formation on its wings is dependent on the second differential of the Ambient temperature around which it is flying. The relation is given by the differential equation: 𝟐

𝒚′′ − 𝒚′ = 𝒕𝟐 𝒆𝒕

Interpret this differential Equation using MATLAB.

MATLAB Code: fun=@(y,t) y2-y1^2=t.^2*exp(t) A= [0 1.41]; B = 5.65; [a,b] = ode23(@(a,b) fun,A, B) plot(a,b,'-o')

Result: a=

1.343 1.2310 1.2230 1.5630

59

0.7800 0.7050 0.8460 0.9870 1.1280 1.2690 1.4100 b= 5.6500 5.6898 5.8090 6.0079 6.2862 6.6441 7.0814 7.5983 8.1948 8.8707 9.6262

60

Graph:

Interpretation: Thus the Deposition of ice gradients for a given range of temperatures are being computed.

3. The visibility of the Aircraft’s Wing mounted lights varies with the change in humidity. The expresssion that gives a relation between the intensity from the light with the humidity is given by: 𝑦′′2 + 3𝑦 = 𝑡 3 − 1 Where y denotes the Visibility in meters and t represents he humidity present in that place. Solve this differential equation using MATLAB.

61

MATLAB Code: fun=@(y,t) y2.^2+3*y=t.^3-1 A= [5 1]; B = 2; [m,n] = ode113(@(x,y) fun,A, B) plot(x,y,'-o')

Result: m= 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 -2.9999 -2.9998 -2.9997 -2.9993

62

-2.9987 -1.9973 -1.9946 -1.9892 -1.9784 -1.9568 -1.9136 -1.8273 -1.6546 -1.3546 -1.0546 -0.7546 -0.4546 -0.1546 0.1454 0.4454 0.7454 1.0000

63

n=

0 -0.0000 -0.0000 -0.0001 -0.0001 -0.0002 -0.0005 -0.0010 -0.0020 -0.0040 -0.0081 -0.0162 -0.0323 -0.0646 -0.1287 -0.2558 -0.5052 -0.9853 -1.8737 -3.3886

64

-5.3233 -6.5786 -7.3495 -7.7827 -7.9760 -7.9787 -7.7918 -7.3672 -6.7500

Interpretation: Hence the variation between these two factors are being computed.

65

Solving Parabolic Partial Differential Equations Theory: A parabolic partial differential equation is a type of partial differential equation (PDE). Parabolic PDEs are used to describe a wide variety of time-dependent phenomena, including heat conduction, particle diffusion, and pricing of derivative investment instruments. To define the simplest kind of parabolic PDE, consider a real-valued function u(x,y) of two independent real variables, x and y. A second-order, linear, constant-coefficient PDE for takes the form

and this PDE is classified as being parabolic if the coefficients satisfy the condition.

Usually x represents one-dimensional position and y represents time, and the PDE is solved subject to prescribed initial and boundary conditions. The name "parabolic" is used because the assumption on the coefficients is the same as the condition for the analytic geometry equation

to define a planar parabola. The basic example of a parabolic PDE is the one-dimensional heat equation,

where u(x,t) is the temperature at time t and at position x along a thin rod, and α is a positive constant (the thermal diffusivity). The symbol Ut signifies the partial derivative of u with respect to the time variable t, and similarly Uxx is the second partial derivative with respect to x. For this example, t plays the role of y in the general second-order linear PDE: A=𝛼, E=1, and the other coefficients are zero.

66

Functions: 1. pdepe : pdepe Solve initial-boundary value problems for systems of parabolic and elliptic PDEs in one space variable and time. In the form of: 𝑐 (𝑥, 𝑡, 𝑢,

𝜕𝑢 𝜕𝑢 𝜕 𝑚 𝜕𝑢 𝜕𝑢 ) = 𝑥 −𝑚 (𝑥 𝑓(𝑥, 𝑡, 𝑢, )) + 𝑠(𝑥, 𝑡, 𝑢, ) 𝜕𝑥 𝜕𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑥

The PDEs hold for t0 ≤ t ≤ tf and a ≤ x ≤ b. The interval [a,b] must be finite. m can be 0, 1, or 2, corresponding to slab, cylindrical, or spherical symmetry, respectively. If m > 0, then a must be ≥ 0. In the Equation, f (x,t,u,∂u/∂x) is a flux term and s (x,t,u,∂u/∂x) is a source term. The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix c (x,t,u,∂u/∂x). The diagonal elements of this matrix are either identically zero or positive. An element that is identically zero corresponds to an elliptic equation and otherwise to a parabolic equation. There must be at least one parabolic equation. An element of c that corresponds to a parabolic equation can vanish at isolated values of x if those values of x are mesh points. Discontinuities in c and/or s due to material interfaces are permitted provided that a mesh point is placed at each interface.

And pdepe syntax are: sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) [sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) Where, m: A parameter corresponding to the symmetry of the problem. m can be slab = 0, cylindrical = 1, or spherical = 2. pdefun: A handle to a function that defines the components of the PDE.

67

icfun : A handle to a function that defines the initial conditions. bcfun: A handle to a function that defines the boundary conditions. xmesh: A vector [x0, x1, ..., xn] specifying the points at which a numerical solution is requested for every value in tspan. The elements of xmesh must satisfy x0 < x1 < ... < xn. The length of xmesh must be >= 3. tspan: A vector [t0, t1, ..., tf] specifying the points at which a solution is requested for every value in xmesh. The elements of tspan must satisfy t0 < t1 < ... < tf. The length of tspan must be >= 3. options: Some options of the underlying ODE solver are available in pdepe: RelTol, AbsTol, NormControl, InitialStep, MaxStep, and Events. In most cases, default values for these options provide satisfactory solutions. See odeset for details. 2. pdeval : pdeval Evaluates numerical solution of PDE using output of pdepe and its syntax is: [uout,duoutdx] = pdeval(m,x,ui,xout) Where , m: Symmetry of the problem: slab = 0, cylindrical = 1, spherical = 2. This is the first input argument used in the call to pdepe. x: A vector [x0, x1, ..., xn] specifying the points at which the elements of ui were computed. This is the same vector with which pdepe was called.

68

ui: A vector sol(j,:,i) that approximates component i of the solution at time tf and mesh points xmesh, where sol is the solution returned by pdepe. xout: A vector of points from the interval [x0,xn] at which the interpolated solution is requested.

Questions: 1. Solve the Heat Equation using MATLAB.

MATLAB Code: mu = 1.0; Umax = 5.0; my = 100; ymax = 10.0; mt = 40; tmax = 100.0; A = tmax/((nt+1)^2); y = linspace(0,ymax,my); t = A*((1:mt)+1).^2; u = pdepe(0,@pdfslpde,@pdfslic,@pdfslbc,y,t,[],mu,Umax);

69

2. Solve the 1-Dimensional Diffusion Equation of Helium gas in Air using MATLAB.

MATLAB Code: qv = 1.0; vmax = 5.0; qy = 100; nmax = 10.0; qt = 40; tmax = 100.0; A = tmax/((qt+1)^2); n = linspace(0,nmax,qn); t = A*((1:qt)+1).^2; v = pdepe(0,@pdfslpde,@pdfslic,@pdfslbc,n,t,[],qu,vmax);