% Script 1 For Week 10 % Create a derivative function of time and place: dydt=@(t,y)f(t,y) % Look at a direction field to see what the solutions will look like. f = @(t,y) -y.^2+t dirfield(f, -1:0.2:4, -2:0.2:2) Euler's Method: y(1) and time vector t given. y(i+1)=y(i)+dydt(t(i),y(i))*(t(i+1)-t(i)); Heun's Method: Use an intermediary estimate to try to improve the step. y(1) and time vector t given. y0(i+1)=y(i)+dydt(t(i),y(i))*(t(i+1)-t(i)); y(i+1)=y(i)+(dydt(t(i),y(i))+dydt(t(i),y0(i)))/2*(t(i+1)-t(i)); https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html Try typing odeexamples on the commandline Various ode solvers: [TOUT,YOUT] = odeXXX(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system of differential equations y' = f(t,y) from time T0 to TFINAL with initial condition(s) Y0 ode45 Solve nonstiff differential equations; medium order method [t,y]=ode45(f,[0 3],[-0.5 0 0.5 1 1.5 2]); dirfield(f, -1:0.2:4, -2:0.2:2), hold on, plot(t,y), hold off % We solve for two different initial conditions for the y variable as a % function of t. y(time 0)=2 and y(time 0)=0. [t,y]=ode45(f,[0 20],[2 0]); figure(3) plot(t,y(:,1)); figure(4) plot(t,y(:,2)); % ode23 Solve non-stiff differential equations, low order method. % vdp1 is the built-in van der Pols equation slope function given by % dydt = [y(2); (1-y(1)^2)*y(2)-y(1)] % ode113 Solve non-stiff differential equations, variable order method. % For stiffer systems (with faster mixed with slower-varying components) % use one of the later ones. Note that they use different numbers of % domain points selected in [0,20]. % You have to really zoom in to distinguish the plots. [t0,y0]=ode45(@vdp1,[0 20],[2 0]); [t1,y1]=ode23(@vdp1,[0 20],[2 0]); [t2,y2]=ode113(@vdp1,[0 20],[2 0]); [t3,y3]=ode15s(@vdp1,[0 20],[2 0]); [t4,y4]=ode23s(@vdp1,[0 20],[2 0]); [t5,y5]=ode23t(@vdp1,[0 20],[2 0]); [t6,y6]=ode23tb(@vdp1,[0 20],[2 0]); % Creates van der Pol solution for mu=1. The first column of each yn % is the solution function, the second is its derivative. The solution starts % out with value 2 and its derivative with value 0. % first few values of t1: % 0 % 0.0000 % 0.0002 % 0.0012 % 0.0062 % 0.0312 % 0.0768 % first few values of y1 = % 2.0000 0 % 2.0000 -0.0001 % 2.0000 -0.0005 % 2.0000 -0.0025 % 2.0000 -0.0124 % 1.9991 -0.0596 % 1.9945 -0.1372 [length(t0),length(t1),length(t2),length(t3),length(t4),length(t5),length(t6)] % Different adaptive methods create different solution lengths % on the same interval. plot(t0,y0,t1,y1,t2, y2,t3, y3,t4, y4, t5, y5,t6, y6) % We now try the same thing with the much stiffer van der Pol with mu=1000. % First create function M-file vanderpol. *************************** yp=vanderpol(t,y,mu) yp=@(y) [y(2);mu*(1-y(1)^2).*y(2)-y(1)]; end **************************** [t0,y0]=ode45(@vanderpol,[0 6000],[1 1],[],1); [t1,y1]=ode23(@vanderpol,[0 6000],[1 1]); [t2,y2]=ode113(@vanderpol,[0 6000],[1 1]); [t3,y3]=ode15s(@vanderpol,[0 6000],[1 1]); [t4,y4]=ode23s(@vanderpol,[0 6000],[1 1]); [t5,y5]=ode23t(@vanderpol,[0 6000],[1 1]); [t6,y6]=ode23tb(@vanderpol,[0 6000],[1 1]); plot(t1,y1(:,1),t2, y2(:,1),t3, y3(:,1),t4, y4(:,1), t5, y5(:,1),t6, y6(:,1)) % Predator-prey model: % Create the function file. Notice there is no time dependence in this % particular DE. The slope fields are NOT time dependent. % The solvers work either way. ************ function yp = predprey(t,y,a,b,c,d) yp = [a*y(1)-b*y(1)*y(2);-c*y(2)+d*y(1)*y(2)]; ************ % and run: h=0.0625;tspan=[0 40];y0=[2 1]; a=1.2;b=0.6;c=0.8;d=0.3; [t y] = eulersys(@predprey,tspan,y0,h,a,b,c,d); subplot(2,2,1);plot(t,y(:,1),t,y(:,2),'--') legend('prey','predator');title('(a) Euler time plot') subplot(2,2,2);plot(y(:,1),y(:,2)) title('(b) Euler phase plane plot') [t y] = rk4sys(@predprey,tspan,y0,h,a,b,c,d); subplot(2,2,3);plot(t,y(:,1),t,y(:,2),'--') title('(c) RK4 time plot') subplot(2,2,4);plot(y(:,1),y(:,2)) title('(d) RK4 phase plane plot') %******************************************************* % The following is the Lorenz Equation example. % This DE is not time dependent either. h=0.03125;tspan=[0 20];y0=[5 5 5]; a=10;b=8/3;r=28; ************* function y = lorenz(t,y,a,b,r) y = [-a*y(1)+a*y(2); r*y(1)-y(2)-y(1)*y(3); -b*y(3)+y(1)*y(2) ]; ************* [t y] = ode45(@lorenz,tspan,y0,h,a,b,r); figure(28) subplot(4,1,4); plot(t,y(:,1),t,y(:,2),t,y(:,3)) title('(d) x, y and z time plot') xlabel('time');ylabel('coordinate values'); legend('x','y','z') legend('location','northwest') subplot(4,1,1);plot(y(:,1),y(:,2)) title('(a) projection: x versus y') xlabel('x');ylabel('y'); subplot(4,1,2);plot(y(:,1),y(:,3)) title('(b) projection: x versus z') xlabel('x');ylabel('z'); subplot(4,1,3);plot(y(:,2),y(:,3)) title('(c) projection: y versus z') xlabel('y');ylabel('z'); figure(29) plot3(y(:,1),y(:,2),y(:,3)) xlabel('x');ylabel('y');zlabel('z');grid title('Lorenz Solution Path') % Now let's take a look at how a tiny modification of % one of the initial conditions affects the solution. y1=[5.001 5 5]; [tchanged ychanged] = ode45(@lorenz,tspan,y1,h,a,b,r); figure(30) plot(t,y(:,1),tchanged,ychanged(:,1)) title('(d) x coordiate versus t') legend('x=y=z=5','x=5.001, y=z=5') legend('location','northwest') % Note that the two solutions coincide for the first 12 seconds or so, then diverge wildly. %******************************************************* function dy = Plinyode(t,y) global siphon Rt = 0.05; r = 0.007; yhi = 0.1; ylo = 0.025; C = 0.6; g = 9.81; Qin = 0.00005; if y(1) <= ylo siphon = 0; elseif y(1) >= yhi siphon = 1; end Qout = siphon * C * sqrt(2 * g * y(1)) * pi * r ^ 2; dy = (Qin - Qout) / (pi * Rt ^ 2); %******************************************************* % This solves the text "Pliny's fountain" problem. tspan = [0 100]; y0 = 0; [tp,yp]=ode45(@Plinyode,tspan,y0); figure(8) plot(tp,yp) title('ode45') ylim([0,0.11]) xlabel('time, (s)') ylabel('water level in tank, (m)') % Let's do the same thing with the ode23s solver. Note big differences! % Who is right? (Probably both wrong in detail for this VERY stiff system.) tspan = [0 100]; y0 = 0; [tp,yp]=ode23s(@Plinyode,tspan,y0); figure(9) plot(tp,yp) title('ode23s') ylim([0,0.11]) xlabel('time, (s)') ylabel('water level in tank, (m)')