numdiff2.m

求解predator-prey的代码

%  numdiff2.m
% Matlab file for Part 2 of the Numerical Solutions
% of Differential Equations module

disp('********************************************')
disp('Part 2:  Numerial Methods in a Helper')
disp('         Application.')
disp('********************************************')
disp('  ')    

    format short

    disp('Step 1: ')
    disp('In order to draw the slope (direction),')
    disp('field for a differential equation,')
    disp('you must first define a MATLAB function')
    disp('giving the right hand side of the')
    disp('differential equation in an m-file')
    disp('called dfun.m.')
    disp(' ')
    disp('------------------------------------------')  
    disp('An example dfun.m file is included')
    disp('with the files provided with this')
    disp('module.  The original is set up for')
    disp('the differential equation')
    disp('   dy/dt = 2 cos(t) - t*y. ')
    disp(' ')
    disp('You can list the dfun.m file while in')
    disp('MATLAB by entering ''''type dfun''''. ')
    disp(' ')
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ')

    keyboard;
    disp(' ')

    disp('------------------------------------------')
    disp('Step 1 (cont.): ')
    disp('The MATLAB function:')
    disp('    slpfield(tmin,tmax,ymin,ymax) ')
    disp('for plotting slope fields is contained')
    disp('in a m-file provided with this module.')
    disp(' ')
    disp('Copy the following command and paste it')
    disp('at a MATLAB prompt, then execute it.')
    disp(' ')
    disp('     slpfield(0,8,-1,3) ')
    disp(' ')
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ')

    keyboard;
    disp(' ')

    disp('-------------------------------------------')
    disp('Step 2: ')
    disp('Set n=20 below and execute the following')
    disp('commands to define equally spaced time values')
    disp('and then use the recursive relation of Euler''s ')
    disp('method to generate y values at these times')
    disp(' ')
    disp('The last commands in the set generate the plot')
    disp('of the (t,y) pairs you have generated and hold')
    disp('the plot so we can add to it. ')
    disp(' ')
    disp('Copy the following commands as a group and paste')
    disp('them at a MATLAB prompt, then execute them.')
    disp(' ')
    disp('  n= ? ')
    disp('  Delta=8/n; ')
    disp('  clear t y ')
    disp('  t=0:Delta:8; ')
    disp('  y0=2; y(1)=y0; ')
    disp('  for k=2:n+1 ')
    disp('    slope(k-1)=dfun(t(k-1),y(k-1));  ')
    disp('    y(k)=y(k-1)+slope(k-1)*Delta; ')
    disp('  end ')
    disp('  plot(t,y,''bo'') ')
    disp('  grid on ')
    disp('  hold on ')
    disp(' ')    
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ')

    keyboard;
    disp(' ')

    disp('-------------------------------------------------')
    disp('Step 3: ')
    disp('Before we add extra graphs to our plot, let''s ')
    disp('save the (t,y) pairs from our last plot so we can ')
    disp('use them later.')
    disp(' ')
    disp('Enter:')
    disp(' ')
    disp(' tE20=t; yE20=y; ')
    disp(' ')
    disp('-------------------------------------------------')
    disp('Now, change the number of steps n to 40 in the')
    disp('commands of step 2, recalculate, and replot.')
    disp('  Change the plot color ''b'' to another color')
    disp('  in the plot command.')
    disp('    After your plot, save the (t,y) pairs as ')
    disp('    tE40 and yE40. ')
    disp(' ')
    disp('Do the same for n=80.  (Don''t forget to change ')
    disp('the plot color again.) ')
    disp(' ')
    disp('------------------------------------------------')
    disp('Describe the changes you see as the number of')
    disp('steps goes from 20 to 40 to 80.')
    disp('Use MATLAB comments in your diary file.')
    disp(' ')
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ') 

    keyboard;
    disp(' ')

    disp('-------------------------------------------------')
    disp('Step 4: ')
    disp('The Improved Euler method uses the Euler Method')
    disp('to predict y(k) at t(k) starting from y(k-1)')
    disp('at t(k-1).  Then the slope is computed at')
    disp('the predicted point (t(k), y(k)).  ')
    disp(' ')
    disp('This new slope at t(k) is averaged with slope(k-1)')
    disp('at t(k-1) and a better estimate of y(k) is')
    disp('computed based on using the average slope')
    disp('to compute the rise. ')
    disp(' ')
    disp('Set n=20 below, then Enter:')
    disp(' ')
    disp('  n= ? ')
    disp('  Delta=8/n; ')
    disp('  clear t y ')
    disp('  t=0:Delta:8; ')
    disp('  y0=2; y(1)=y0; ')
    disp('  for k=2:n+1 ')
    disp('    slope1(k-1)=dfun(t(k-1),y(k-1));  ')
    disp('    y1(k)=y(k-1) + slope1(k-1)*Delta; ')
    disp('    slope2(k)=dfun(t(k),y1(k)); ')
    disp('    aveslope=0.5*( slope1(k-1)+slope2(k) );')
    disp('    y(k)=y(k-1) + aveslope*Delta; ')
    disp('  end ')
    disp('  plot(t,y,''bo'') ')
    disp('  grid on ')
    disp('  hold on ')
    disp(' ')
    close    
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ')
    

    keyboard;
    disp(' ')

    disp('-------------------------------------------------')
    disp('Step 4 (cont.): ')
    disp('Before we add extra graphs to our plot, let''s ')
    disp('save the (t,y) pairs from our last plot so we can ')
    disp('use them later.  Use ''''IE'''' for Improved Euler.')
    disp(' ')
    disp('Enter:')
    disp(' ')
    disp(' tIE20=t; yIE20=y; ')
    disp(' ')
    disp('-------------------------------------------------')
    disp('Now, change the number of steps n to 40 in the')
    disp('commands of step 4, recalculate, and replot.')
    disp('  Change the plot color ''b'' to another color')
    disp('  in the plot command.')
    disp('    After your plot, save the (t,y) pairs as ')
    disp('    tIE40 and yIE40. ')
    disp(' ')
    disp('Do the same for n=80.  (Don''t forget to change ')
    disp('the plot color again and save the (t,y) values.) ')
    disp(' ')
    disp('------------------------------------------------')
    disp('Describe the changes you see as the number of')
    disp('steps goes from 20 to 40 to 80.')
    disp('Use MATLAB comments in your diary file.')
    disp(' ')
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ')

    keyboard;
    disp(' ')
 
    disp('--------------------------------------------')
    disp('Step 5: ')
    disp('MATLAB has a built in numerical DE solver called')
    disp('ode45 that uses a 4th order Runge-Kutta Method.')
    disp('Let''s show the RK solution on top of the ')
    disp('direction field. ')
    disp(' ')
    disp('Enter:')
    disp(' ')
    disp('  slpfield(0,8,-1,3); hold on ')
    disp('  [tRK,yRK]=ode45(''dfun'', [0,8], y0); ')
    disp('  plot(tRK,yRK,''k'') ')
    disp(' ')
    close
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ')

    keyboard;
    disp(' ')

    disp('--------------------------------------------')
    disp('Step 5 (cont.) ')
    disp('Since we have saved all the points in our plots')
    disp('we can easily compare them.')
    disp(' ')
    disp('To compare the n=20 Euler plot with the Runge-Kutta')
    disp('''''exact'''' plot.....Enter: ')
    disp(' ')
    disp('  clf   % Clears figure ')
    disp('  plot(tRK,yRK,''k''); grid on; hold on ')
    disp('  plot(tE20,yE20,''bo'') ')
    disp(' ')
    disp('-----------------------------------------------')
    disp('All you have to do to compare other plots to ')
    disp('the RK plot, is change tE20,yE20 to other ')
    disp('Euler or Improved Euler sets of (t,y) pairs. ')
    disp(' ')
    disp('-----------------------------------------------')
    disp('Use MATLAB comments to record your observations')
    disp('and conclusions in your diary file.')
    disp(' ')
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ')

    keyboard;
    disp(' ')

    disp('--------------------------------------------')
    disp('Step 6:')
    disp('Use the commands in step 2 to calculate the')
    disp('Euler plot with n=10 steps.   ')
    disp(' ')
    disp('--------------------------------------------')
    disp('Describe what is happening and why. ')
    disp('Use MATLAB comments in your diary file. ')
    disp(' ')
    disp(' ')
    disp('When you have finished part 2, to go on')
    disp('to part 3 of this module, type: numdiff3')
    disp(' ')