1 / 7

Lecture 29: Numerical integration of Ordinary Differential equations

Lecture 29: Numerical integration of Ordinary Differential equations. Euler’s method. euler1a.m. euler1atest.m. %Euler's Method for Solving Initial Value Problems %Use with ydot.m to evaluate rhs of differential equation % Input: interval [a,b], initial value y0, step size h

dixon
Download Presentation

Lecture 29: Numerical integration of Ordinary Differential equations

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. Lecture 29: Numerical integration of Ordinary Differential equations Euler’s method. euler1a.m euler1atest.m

  2. %Euler's Method for Solving Initial Value Problems %Use with ydot.m to evaluate rhs of differential equation % Input: interval [a,b], initial value y0, step size h % Output: time steps t, solution y % Example usage: [t y]=euler1a(func,[0 1],1,0.1); function [t y]=euler1a(func,int,y0,h) clear t; clear y; t(1)=int(1); y(1)=y0; n=round((int(2)-int(1))/h); for i=1:n t(i+1)=t(i)+h; y(i+1)=y(i)+h*func(t(i),y(i)); end

  3. %euler1atest.m - test of function euler1a. Runing require file euler1a.m % Example usage: y=euler1a(func,[0 1],1,0.1); f1=inline('t*y','t','y') %definition of inline function with 2 arguments t and y y1exact=inline('exp(t.^2/2)'); [t y]=euler1a(f1,[0 1], 1, 0.1); y1exact(t); [tode45 yode45]=ode45(f1,[0 1], 1); plot(t,y,'r',tode45,yode45,'b'); legend('euler','exact'); ii=1:10; hvec=0.1./2.^ii; for i=1:length(hvec) [t y]=euler1a(f1,[0 1], 1, hvec(i)); disp(['h=',num2str(hvec(i)),' Error in eulers method =',num2str(y(length(y))-yode45(length(yode45)))]); %display value of error end

  4. >> euler1atest f1 = Inline function: f1(t,y) = t*y h=0.05 Error in eulers method =-0.05278 h=0.025 Error in eulers method =-0.026921 h=0.0125 Error in eulers method =-0.013598 h=0.00625 Error in eulers method =-0.0068341 h=0.003125 Error in eulers method =-0.0034259 h=0.0015625 Error in eulers method =-0.0017152 h=0.00078125 Error in eulers method =-0.00085815 h=0.00039063 Error in eulers method =-0.00042921 h=0.00019531 Error in eulers method =-0.00021464 h=9.7656e-005 Error in eulers method =-0.00010733

  5. Euler’s method for systems euler2.m % Program 6.2 Vector version of Euler method % Inputs: interval [a,b], initial vector y0, step size h % Output: time steps t, solution y % Example usage: y=euler2([0 1],[0 1],0.1); function [t,y]=euler2(int,y0,h) t(1)=int(1); y(1,:)=y0; n=round((int(2)-int(1))/h); for i=1:n t(i+1)=t(i)+h; y(i+1,:)=eulerstep(t(i),y(i,:),h); end plot(t,y(:,1),t,y(:,2)); function y=eulerstep(t,y,h) %one step of the Euler method %Input: current time t, current vector y, step size h %Output: the approximate solution vector at time t+h y=y+h*ydot(t,y); function z=ydot(t,y) z(1) = y(2)^2-2*y(1); z(2) = y(1)-y(2)-t*y(2)^2;

  6. >> euler2([0 1],[0 1],0.1) t = Columns 1 through 10 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 Column 11 1.0000 y = 0 1.0000 0.1000 0.9000 0.1610 0.8119 0.1947 0.7336 0.2096 0.6636 0.2117 0.6006 0.2054 0.5437 0.1939 0.4921 0.1793 0.4453 0.1633 0.4029 0.1469 0.3643

  7. Euler’s method for system of 4 equations orbit.m %Program 6.4 Plotting program for one-body problem % Inputs: int=[a b] time interval, initial conditions % ic = [x0 vx0 y0 vy0], x position, x velocity, y pos, y vel, % h = stepsize, p = steps per point plotted % Calls a one-step method such as trapstep.m % Example usage: orbit([0 100],[0 1 2 0],0.01,5) function z=orbit(int,ic,h,p) n=round((int(2)-int(1))/(p*h)); % plot n points x0=ic(1);vx0=ic(2);y0=ic(3);vy0=ic(4); % grab initial conds y(1,:)=[x0 vx0 y0 vy0];t(1)=int(1); % build y vector set(gca,'XLim',[-5 5],'YLim',[-5 5],'XTick',[-5 0 5],'YTick',... [-5 0 5],'Drawmode','fast','Visible','on','NextPlot','add'); cla; sun=line('color','y','Marker','.','markersize',25,... 'xdata',0,'ydata',0); drawnow; head=line('color','r','Marker','.','markersize',25,... 'erase','xor','xdata',[],'ydata',[]); tail=line('color','b','LineStyle','-','erase','none',... 'xdata',[],'ydata',[]); %[px,py,button]=ginput(1); % include these three lines %[px1,py1,button]=ginput(1); % to enable mouse support %y(1,:)=[px px1-px py py1-py]; % 2 clicks set direction for k=1:n for i=1:p t(i+1)=t(i)+h; y(i+1,:)=eulerstep(t(i),y(i,:),h); end y(1,:)=y(p+1,:);t(1)=t(p+1); set(head,'xdata',y(1,1),'ydata',y(1,3)) set(tail,'xdata',y(2:p,1),'ydata',y(2:p,3)) drawnow; End function y=eulerstep(t,x,h) %one step of the Euler method …. >> orbit([0 100], [0 1 2 0],0.01,5) Show animation

More Related