10 34 Fall 2006 Homework 1 Due Date Wednesday Sept 13th 2006 9 AM Problem 1 Bessel Functions Bessel functions are commonly encountered in heat and mass transfer problem where the geometry is cylindrical Bessel functions of the first kind J x and second kind Y x are the two linearly independent solutions of the differential equation x2 d2y dy x x2 2 y 0 2 dx dx The solution of the equation is exactly J x for boundary condition y 0 1 and dy 0 0 dx This second order differential equation can be converted into a system of two coupled first order differential equations by defining new variable u1 and u2 as follows u1 y u2 dy dx The two differential equations thus obtained are du1 u2 dx du2 u 2 2 1 2 u1 dx x x 1 A matlab code is presented below which solves the above problem with boundary conditions u1 0 1 and u2 0 0 du2 in equation 1 becomes singular To write the matlab code we dx have used the property of Bessel functions J 0 0 5 for 1 Notice that at x 0 Cite as William Green Jr course materials for 10 34 Numerical Methods Applied to Chemical Engineering Fall 2006 MIT OpenCourseWare http ocw mit edu Massachusetts Institute of Technology Downloaded on DD Month YYYY problem 1 HW set 1 Solution of bessel s equation function using ode23s function plot bessel using ode x end u init is the intial condition at t 0 u init 1 0 solve the differential equation using ode23s to generate vectors for x and u x u ode23s diff 0 x end u init plot x u 1 title Bessel function of first kind xlabel x ylabel y x return function diff evaluates the derivatives of u1 ad u2 function f diff x u f zeros size u f 1 u 2 if x 0 f 2 0 5 else f 2 u 2 x u 1 end return 1 Copy the matlab code and run it to generate the plot of Bessel functions Hint Make sure that the file is saved in file name plot bessel usinig ode m although this is not essential it is a good practice to name file names the same as the function they contain Also make sure that the matlab present directory is the same as the directory which contains the file plot bessel using ode m Type the command plot bessel using ode on the matlab prompt 2 The next task is to plot the function J 0 x using the built in matlab command besselj First make a vector of x using the command linspace to generate 40 equally spaced values between 0 and 10 Then use a for loop to iterate over all Cite as William Green Jr course materials for 10 34 Numerical Methods Applied to Chemical Engineering Fall 2006 MIT OpenCourseWare http ocw mit edu Massachusetts Institute of Technology Downloaded on DD Month YYYY the values of x to calculate J 0 x using the command besselj Look at matlab help to figure out how to use besselj and linspace commands 3 If you looked at the help for besselj you will realize that value of J 0 x can be calculated using the infinite series 1 k x 2 2 k k 0 k k 1 J x The good news for Bessel functions is that each consecutive term of the Bessel Series quickly approaches 0 Write a fuction my bessel nu x which takes in a value of and x respectively and returns the value of J x and also the number of terms in the series used to calculate it Then plot the value of J x for x 0 10 and compare the graph to the ones that you got in part 1 and part 2 Also plot the number of terms n required to achieve a converged value of J x against x HINT Let 1 k x 2 2 k tk then k k 1 J x tk but for our purposes we will approximate J x with the k 0 expression n J x tk k 0 where n is sufficiently big We will say that the value of J x has converged when tn 1 n t k 0 0 001 k You might want to use the while loop in matlab to check when the n 1 th term becomes less than 0 1 of the overall sum Cite as William Green Jr course materials for 10 34 Numerical Methods Applied to Chemical Engineering Fall 2006 MIT OpenCourseWare http ocw mit edu Massachusetts Institute of Technology Downloaded on DD Month YYYY Problem 2 Equation Solvers Cubic equations of state are often good at describing vapor and liquid behavior but it is typically difficult to find the volume roots analytically Use the fzero fsolve and fminsearch Matlab functions to determine the vapor volume root of the following equation of state for water Redlich Kwong at P 1 atm and T 300 K P a T RT V b V b V b with a T Tr RTc 2 and Pc b RTc Pc where Tr Tr 1 2 1 0 0 42748 0 08664 for water Tc 647 1 K Pc 217 7 atm Compare the performance of these solvers over a range of initial volume guesses from 1 to 50 L mole Plot the values returned for each solver as a function of the initial guess Do some solvers appear to be more robust than others Some useful Matlab tools for this problem fzero fsolve fminsearch global linspace plot function Cite as William Green Jr course materials for 10 34 Numerical Methods Applied to Chemical Engineering Fall 2006 MIT OpenCourseWare http ocw mit edu Massachusetts Institute of Technology Downloaded on DD Month YYYY Problem 3 Reading and Writing Files In this problem you will have to use some of the I O utilities in Matlab that allow one to read and write text and Excel files Write a Matlab script that reads the text and excel files posted on the website p2 text txt and p2 excel xls containing the modified Arrhenius parameters of a reaction and plots y axis log scale the rate constant over the temperature range from 300 K to 1500 K Also print the rate constant in intervals of 100 K to both a text file and an Excel file in the following form Rate constant data for NH2OH HO2 NH2O HOOH Temp K 300 1500 k T cm3 mole s XXX YYY Some useful Matlab tools for this problem textscan xlsread fprintf xlswrite char fopen fclose length strmatch for plot Cite as William Green Jr course materials for 10 34 Numerical Methods Applied to Chemical Engineering Fall 2006 MIT OpenCourseWare http ocw mit edu Massachusetts Institute of Technology Downloaded on DD Month YYYY