numdiff3.m

求解predator-prey的代码

%  numdiff3.m
% Matlab file for Part 3 of the Numerical Solutions
% of Differential Equations module

disp('********************************************')
disp('Part 2:  The Logistic Equation ')
disp('********************************************')
disp('  ')    

    format short

    disp('Step 1: ')
    disp('In order to study the logistic equation,')
    disp('you first need to use a text editor to')
    disp('change the right hand side of the DE')
    disp('given in the file dfun.m provided with')
    disp('this module.')
    disp(' ')
    disp('Edit the file dfun.m so the the new DE is')
    disp(' ')
    disp('    dy/dt = y.*(1-y)  ')
    disp(' ')
    disp('  (Don''t forget the period after the')
    disp('   first y on the right hand side.')
    disp(' ')
    disp('------------------------------------------')  
    disp('After you have edited the file, you can enter')
    disp('''''type dfun'''' for a check. ')
    disp(' ')
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ')

    keyboard;
    disp(' ')

    disp('------------------------------------------')
    disp('Step 1 (cont.): ')
    disp('We will solve the Logistic differential equation ')
    disp('numerically using the Runge-Kutta solver ode45 ')
    disp('and plot the solution on the direction field.')
    disp(' ')
    disp('Enter: ')
    disp(' ')
    disp('  clf   % Clears figure ')
    disp('  slpfield(0,8,0,1.5); hold on ')
    disp('  y0=0.1  ')
    disp('  [tRK,yRK]=ode45(''dfun'', [0,8], y0); ')
    disp('  plot(tRK,yRK,''r-o'') ')
    disp(' ')
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ')

    keyboard;
    disp(' ')

    disp('-------------------------------------------')
    disp('Step 2: ')
    disp('Let''s use the symbolic capabilities of MATLAB')
    disp('to solve the initial value problem:')
    disp('  dy/dt = y*(1-y), y(0), y''(0)=0 ')
    disp('symbolically. ')
    disp(' ')
    disp('Enter: ')
    disp(' ')
    disp('  syms t ')
    disp('  ysol=dsolve(''D1y = y*(1-y)'',''y(0)=y0'') ')
    disp('  ysol=simple(ysol)')
    disp('  pretty(ysol)  % Nicer format ')
    disp(' ')
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ')

    keyboard;
    disp(' ')

    disp('-------------------------------------------------')
    disp('Step 2 (cont.): ')
    disp('Now, plug the values of the initial condition into')
    disp('the solution and plot it on top of our previously ')
    disp('held plot of the numerical solution. ')
    disp('Enter: ')
    disp(' ')
    disp(' y0 ')
    disp(' y=subs(ysol)')
    disp(' ezplot(y,[0,8]) ')
    disp(' axis([0,8, 0,1.5)] ')
    disp(' ')
    disp('--------------------------------------------')
    disp('Does your symbolic solution match the numerical')
    disp('solution?  Discuss in your diary file.')
    disp(' ')
    disp('To continue afterwards, type the word return')
    disp('and hit enter!')
    disp(' ')
    
    keyboard;
    disp(' ')
    
    disp('--------------------------------------------')
    disp('Step 3: ')
    disp('Repeat steps 1 and 2 with the new initial ')
    disp('condition y0=1.5.')
    disp(' ')
    disp('--------------------------------------------')
    disp('Does your symbolic solution match the numerical')
    disp('solution?  Discuss in your diary file.')
    disp(' ')
    disp(' ')
    disp('When you have finished part 3, to go on')
    disp('to part 4 of this module, type: numdiff4')
    disp(' ')