1 / 11

How to solve ODE's in Matlab

How to solve ODE's in Matlab. Quite easy, once you get the hang of it!. Matlab has a series of solvers. Most often, it is OK to us ODE23 or ODE45 ode23 better for less smooth solutions ode45 better for smoother solutions

lamya
Download Presentation

How to solve ODE's in Matlab

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. How to solve ODE's in Matlab Quite easy, once you get the hang of it!

  2. Matlab has a series of solvers. • Most often, it is OK to us ODE23 or ODE45 • ode23 better for less smooth solutions • ode45 better for smoother solutions • "doc ode23" gives you information on all of them, but requires more math than I have presented to fully understand

  3. ode45 and its friends are called as follows: • [T,Y] = ode45(@odefun,tspan,y0) • odefun(t,y) is a function! • Inputs: • t is the time, a scaler (explain scaler) • y is a column vector of the variables you are solving for • Output: • dy_dt, a column vector of the time derivative of y

  4. Prey are born Prey are eaten, and turned into Carnivores Carnivores die of old age • Example, from our old prey/carnivore system: • y=[P;C] (why the ";"?) • dy_dt is

  5. So what does the code for the derivative function look like? function dPopdt=dPop_dt(t,Pop); %calculate the Lotka-Volterra carnivore prey equations % % d(Prey)/dt = Prey*(alpha-beta*Carnivore) % d(Carnivore)/dt = -Carnivore*(gamma-delta*Prey) % % where % % The state variable is Pop(1) is prey, Pop(2) is carnivore gamma=0.05; delta=0.001; alpha=0.04; beta=0.00075; dPopdt=[1;1]; %Explain this line! dPopdt(1)=Pop(1)*(alpha-beta*Pop(2)); dPopdt(2)=-Pop(2)*(gamma-delta*Pop(1));

  6. And we would call this with: • [T,Pop] = ode45(@dPop_dt,tspan,y0) • tspan=[tstart,tmin] • y0 is a column vector of the initial condition • ode45 returns • T, all the time it computed a solution for • Y, a matrix that has the same number of columns as y0 has rows and the same number of rows as T, which is the solution • each row corresponds to a time in T • each column corresponds to a variable • Yes, this is stupid. (why do I say this?)

  7. Example: P0=[10;10]; %intial condition %use ode45 to solve [tvec,Pvec]=ode45(@dPop_dt,[0,500],P0); plot(tvec,Pvec(:,2),'r*-',tvec,Pvec(:,1),'g-','linewidth',2) xlabel('time, days') ylabel('critter numbers') title('Red are Predators, Green are Prey') grid on

  8. In the lab we will study a radioactive decay series, were one element decays into another which decays into another: • What are the equations for this, given • Element E1, a fraction  of which decays each year to E2, • Element E2, a fraction  of which decays each year to E3, • and etc. • In units of moles (why?)

  9. At “equilibrium”, what would the ratio of E1 to E2 be, roughly? • -- Why?

  10. How to draw a legend on your plots • han=plot(t,x1,’r-’,t,x2,’b-’) • legend(han,’the red line’,’the blue line’)

More Related