80 likes | 276 Views
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
E N D
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 % 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
%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
>> 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
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;
>> 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
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