VIEWS: 35 PAGES: 5 POSTED ON: 3/2/2010 Public Domain
Plotting A Direction Field (Grain Plot) Using MATLAB 1. Write the first-order ODE in the form dx f (t , x) dt Here x(t ) is the unknown function or dependent variable, and t is its argument, or the independent variable. 2. Write a function file using the independent variable, t , and the dependent variable, x , as inputs, and giving the function values f (t , x) as output. Call it, for example, grainfun.m (or any other name you think is appropriate). Example: For the ODE dx x exp(t ) dt you would write the following function file, called grainfun.m: function f = grainfun(t, x) % f = grainfun(t, x) % function to produce direction field for f = dx/dt % using given t and x values f = x*exp(t) return 3. Use the pre-written functions odegr.m (which uses mode23.m) to produce a direction field plot. Example: For the ODE given above, we can produce a direction field plot as follows (assuming you have already written the function file grainfun.m) >> odegr(‘grainfun’, 0, 1, -6, 6) Implementing Euler’s Method To Solve A First-Order IVP 1. Write the IVP in the form dx f (t , x) x(0) 1 dt 2. Write a function file using the independent variable, t , and the dependent variable, x , as inputs, and giving the function values f (t , x) as output. Call it, for example, rhs_function.m (or any other name you think is appropriate). Example: For the ODE dx x exp(t ) dt you would write the following function file, called rhs_function.m: function f = rhs_function(t, x) % f = rhs_function(t, x) % returns the rhs of the ODE for given t and x % values f = x*exp(t) return 3. Use the pre-written function file euler.m to compute a numerical approximation to a set of solution values, and plot them to produce an (approximate) solution curve. Example: For the ODE given above, we can implement Euler’s method as follows (assuming you have already written the function file rhs_function.m) >> [x, t] = euler(10, 0.1, 1) Listing of M-File for Implementing the Euler Method function [x,t] = euler(n,dt,x0) % [x,t] = euler(n,dt) % returns the first n+1 euler approximations using time % step dt and initial condition x0, with initial time % zero. x = x0; t = 0; for i = 2:n+1 % What is the next line doing? % t(i) = (i-1)*dt; % What is the next line doing? % F = rhs_function(t(i-1),x(i-1)); % What is the next line doing? % x(i) = x(i-1) + dt*F; end plot(t,x) return Table of Values Independent Variable Dependent Variable t0 0 x0 t1 t0 t x1 x0 f (t0 , x0 )t t t2 t1 t x2 x1 f (t1 , x1 )t 2t t3 t 2 t x3 x2 f (t2 , x2 )t 3t tn tn 1 t xn xn 1 f (tn 1 , xn 1 )t n t