Problem 2 – Equation Solvers Problem 3 – Reading and Writing Files10.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 ( ) and second kind ( ) are the two linearly independent solutions of the differential equation ( )Jxν()Yxν2222()dy dyxxxydx dxν++− =20 The solution of the equation is exactly for boundary condition and ( )Jxν(0) 1y =(0) 0dydx=. 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: 1uy= 2dyudx= The two differential equations thus obtained are 12duudx= 222121du uudx x xν⎛⎞=− − −⎜⎟⎝⎠ (1) A matlab code is presented below which solves the above problem with boundary conditions. 1(0) 1u = and 2(0) 0u = (Notice that at x =0, 2dudx in equation 1 becomes singular. To write the matlab code we have used the property of Bessel functions ' (0) 0.5Jν=, for 1ν≥ ). 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( (x)'); 'yreturn; %/////////////////////////////////////////////////////////////// %/////////////////////////////////////////////////////////////// %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 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 0()JxCite 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 using the command besselj. [Look at matlab help to figure out how to use besselj and linspace commands]. 0()Jx 3. If you looked at the help for besselj you will realize that value of can be calculated using the infinite series: 0()Jx 20(1)(/2)()!( 1)kkkxJxkkννν+∞=−=Γ++∑ 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 and also the number of terms in the series used to calculate it. Then plot the value of 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 against x. ( )Jxν( )Jxν( )Jxν [HINT: Let, 2(1)(/2)!( 1)kkkxtkkνν+−=Γ++, then 0()kkJxν∞==∑t, but for our purposes we will approximate with the expression, ( )Jxν 0()nkkJx tν=≈∑, where n is sufficiently big. We will say that the value of has converged when: ( )Jxν 100.001nnkktt+=≤∑. 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. ()()( )()()( )()122::1 0 0.42748 0.08664:647.1 217.7rccccrrccaTRTPVb V bV bwithTRTRTa T and bPPwhereTTfor waterTKPatmεσαψασεψ−=−−++====== Ω=== 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) k(T) (cm3/mole-s) 300 XXX … … 1500 YYY Some useful Matlab tools for this problem: textscan, xlsread, fprintf, xlswrite, char, fopen, fclose, length, strmatch, for, plot Cite