% diary for nov30 % solving ODEs % look for ODE solvers in helpwin under funfun >> helpwin % Note that all of them are for solving a first order ODE initial value % problem, y'(t) = f(t,y(t)) on a <= t <= b, y(a) = ya , for some % function f = f(t,z) but do allow y = y(t) to be a vector. % Example: y'' = -y on [0 .. pi], y(0) = 0, y'(0) = 1. % Convert to first order ODE: Z := [y;y'], hence >> f = inline('[z(2); -z(1)]','t','z') f = Inline function: f(t,z) = [z(2); -z(1)] >> a = 0; b = pi; [t,Z] = ode45(f,[a b],[0;1]); plot(t,Z) >> size(Z) ans = 61 2 % Note that the initial value, Z(0) = [0;1] , is required to be put in as % a 1-column matrix, yet the output, xsol, has Z(i,:) as the value of the % solution at t(i), all i. Presumably, matlab made this inconsistent choice % because that makes the plot command work fine, plotting both components of % the solution in one plot statement. % student's question: what would happen if we specified both the initial % value and the function f itself as being row-vector valued?? % Let's try it: >> f = inline('[z(2) -z(1)]','t','z') f = Inline function: f(t,z) = [z(2) -z(1)] >> ode45(f,[a b],[0 1]); ??? Error using ==> upper Function 'upper' not defined for variables of class 'inline'. Error in ==> C:\MATLABR11\toolbox\matlab\funfun\ode45.m On line 290 ==> error([upper(odefile) ' must return a column vector.']) % ... We get back an incomprehensible error message. In short: don't do it. % Which of the two curves is y ?? >> hold on, plot(t,Z(:,1),'linew',2) % The plot confirms that the solution is Z(t) = [sin(t); cos(t)], i.e., % y = sin . % % next problem is a scalar first-order problem: >> f = 'sqrt(t^2 + z^2)-1'; a = -1; b = 4; >> ya = -2; % we start off by plotting the direction field. For that, we need to % evaluate the slope function at many points: >> ftz = inline(f,'t','z') ftz = Inline function: ftz(t,z) = sqrt(t^2 + z^2)-1 >> [tt, zz] = ndgrid(a:.25:b, -2:.25:2); % better vectorize the function: >> vftz = vectorize(ftz) vftz = Inline function: vftz(t,z) = sqrt(t.^2 + z.^2)-1 % ... >> slopes = vftz(tt,zz); >> clf >> dirfield(tt,zz,slopes,.05) % here, for the record, is the m-file dirfield, along with the m-file segplot % it makes use of: >> type dirfield function dirfield(x,y,slopes,s) %DIRFIELD plot direction field % % dirfield(x,y,slopes,s) % % plots straight-line segment of length 2*S, centered at (X(i),Y(i)) and % of slope SLOPES(i), all i. If the points are to be on a grid, put them % there using NDGRID. % % See also SEGPLOT % cb 1dec96; 16mar97; 17nov98 n = prod(size(x)); x = reshape(x,1,n); y = reshape(y,1,n); slopes = reshape(slopes,1,n); direct = [ones(1,n);slopes]; sqrt(1+slopes.^2); direct = s*direct./ans([1 1],:); segplot([x-direct(1,:);y-direct(2,:);x+direct(1,:);y+direct(2,:)]); >> type segplot function segsout = segplot(s,arg2,arg3) %SEGPLOT a collection of segments (to be) plotted as ONE curve. % % segsout = segplot(s,arg2,arg3) % % returns the appropriate sequences SEGSOUT(1,:) and SEGSOUT(2,:) % (containing the segment endpoints properly interspersed with NaN's) % so that PLOT(SEGSOUT(1,:),SEGSOUT(2,:)) plots the straight-line % segment(s) with endpoints (S(1,:),S(2,:)) and (S(d+1,:),S(d+2,:)) , % with S of size [2*d,:]. % % If there is no output argument, the segment(s) will be plotted in % the current figure (and nothing will be returned), % using the linestyle and linewidth optionally specified by ARG2 and ARG3 % as a string and a number, respectively. % cb dec96, mar97, sep97; DIRFIELD is a good test for it. [twod,n] = size(s); d = twod/2; if d<2, error('the input must be of size (2d)-by-n with d>1.'), end if d>2, s = s([1 2 d+1 d+2],:); end NaN; [s; ans(1,ones(1,n))]; segs = [reshape(ans([1 3 5],:),1,3*n); reshape(ans([2 4 5],:),1,3*n)]; if nargout>1 segsout = segs; else symbol=[]; linewidth=[]; for j=2:nargin eval(['arg = arg',num2str(j),';']) if isstr(arg), symbol=arg; else [ignore,d]=size(arg); if ignore~=1, error(['arg',num2str(j),' is incorrect.']) else linewidth = arg; end end end if isempty(symbol) symbol='-'; end if isempty(linewidth) linewidth=1; end plot(segs(1,:),segs(2,:),symbol,'linewidth',linewidth) % plot(s([1 3],:), s([2 4],:),symbol,'linewidth',linewidth) % would also work, without all that extra NaN, but it's noticeably slower. end %%% back to our calculation. We now plot the solution to our ODE on top % of the dirfield: >> hold on % .... for this, we fix the current scaling of the plot, by the command >> axis(axis) >> [tsol,ysol] = ode45(ftz,[a b],ya); plot(tsol,ysol,'r','linew',2) % Now we construct a Taylor series solution, % y(t) = y(a) + y'(a)(t-a) + (y''(a)/2)(t-a)^2 + (y'''(a)/6)(t-a)^3 + ... % term by term, using the differential equation itself to provide us with % the needed derivatives of y at a: % We start with the linear Taylor polynomial >> taylor = ya + ftz(a,ya)*(tsol-a); plot(tsol,taylor) %% Well, we get a straight line that is tangent to the solution at t = a, % but then quickly diverges from it. We hope to do better by adding terms: %% To get y''(a), we differentiate y'(t) = f(t,y(t)), but prefer to do this % by machine, making use of matlab's symbolic toolbox command diff : >> diff(f,'t') ans = 1/(t^2+z^2)^(1/2)*t >> diff(f,'z') ans = 1/(t^2+z^2)^(1/2)*z %% Here is an m-file that starts off with a string for f(t,z) and a string % zprime for dz/dt and, from it, fashions the string that gives % d f(t,z(t))/dt, hence can be made into function using the inline command` >> type dtftz function deriv = dtftz(f,zprime) %DTFTZ total t derivative of f(t,z) for given dz/dt % % deriv = dtftz(f,zprime) % % returns the string containing the total derivative % % deriv = D_t f(t,z(t)) + z' D_z f(t,z) % % wrto t of the function f(t,z) described by the string f , % using z'(t) as specified by zprime . deriv = char(simplify(sym( [ ... char(diff(f,'t')), '+(', zprime, ')*(', char(diff(f,'z')), ')'... ]))); %%%%% note that the output from diff is a symbolic expression rather than a %% string: >> diff(f,'t') ans = 1/(t^2+z^2)^(1/2)*t >> whos ans Name Size Bytes Class ans 1x1 162 sym object Grand total is 20 elements using 162 bytes % ... hence for our purposes, we need to convert it into string again, using % the command char. With that, we can concatenate the strings for diff(f,'t'), for zprime multiplying diff(f,'z'), and, in general, that will % be a horrendous expression. So, we use the symbolic toolbox again to % simplify this expression, -- except that the command simplify expect a % symbolic expression rather than a string, so you see conversion to symbolic,k % followed by simplify, followed by char, to make the final output a string % again. % Ready for use?? >> f2 = dtftz(ftz,f) ??? Error using ==> diff Function 'diff' not defined for variables of class 'inline'. Error in ==> c:\courses\412\00\dtftz.m On line 13 ==> deriv = char(simplify(sym( [ ... %% Well, I forgot that dtftz expects strings, not inline functions! >> f2 = dtftz(f,f) f2 = (t+z*(t^2+z^2)^(1/2)-z)/(t^2+z^2)^(1/2) % ... but, for evaluation, we need to convert to inline function: >> f2tz = inline(f2,'t','z'); % We use it just at a, to get y''(a) = f2(a,y(a)) >> taylor = taylor + f2tz(a,ya)/2*(tsol-a).^2; plot(tsol,taylor,'k') %% now our approximate solution stick with the solution a little while longer % but then diverges in the opposite direction. % We try for the next term: >> f3 = dtftz(f2,f) f3 = (t^2+z^2-(t^2+z^2)^(1/2)*z*t+2*z*t+t^4+2*t^2*z^2+z^4-2*(t^2+z^2)^(1/2)*t^2-(t^2+z^2)^(1/2)*z^2)/(t^2+z^2)^(3/2) %% (wouldn't want to have done this by hand :-) >> f3tz = inline(f3,'t','z'); >> taylor = taylor + f3tz(a,ya)/6*(tsol-a).^3; plot(tsol,taylor,'g') %% Ok, better yet, but only for a little bit longer. %% MORAL: polynomials are only good for LOCAL approximation. For GLOBAL % approximation, use PIECEWISE polynomials! %% One-step methods: % Here is an m-file for experimenting with low-order one-step methods: >> type show_ode function [tout,yout] = show_ode(f,a2b,ya,n,method,color) %SHOW_ODE sample onestep methods % % [t,y] = show_ode(f,[a,b],ya,n,method) a = a2b(1); b = a2b(2); t = linspace(a,b,n+1); h = (b-a)/n; y = zeros(length(ya),length(t)); y(:,1) = ya; for j=1:n switch method(1:3) case 'eul' dely = h*f(t(j),y(j)); case 'rk2' F1 = h*f(t(j),y(j)); F2 = h*f(t(j+1), y(:,j) + F1); dely = (F1 + F2)/2; case 'rk4' F1 = h*f(t(j),y(j)); F2 = h*f(t(j)+h/2, y(:,j) + F1/2); F3 = h*f(t(j)+h/2, y(:,j) + F2/2); F4 = h*f(t(j+1),y(:,j)+F3); dely = (F1 + 2*(F2+F3)+F4)/6; otherwise, error(['specified method: ''', method, ''' not recognized.']) end y(:,j+1) = y(:,j) + dely; end if nargout==0 if nargin<6, color = 'g'; end plot(t,y,color,'linew',2) else tout = t; yout = y; end >> show_ode(f,[a b],ya,10,'euler') Warning: Subscript indices must be integer values. > In c:\courses\412\00\show_ode.m at line 13 ??? Index into matrix is negative or zero. See release notes on changes to logical indices. Error in ==> c:\courses\412\00\show_ode.m On line 13 ==> dely = h*f(t(j),y(j)); %% what happened?? A student kindly pointed out that ode-solvers, including % this one, expect to be given the slope function f as a FUNCTION, not a % string. >> show_ode(ftz,[a b],ya,10,'euler') >> show_ode(ftz,[a b],ya,25,'euler') >> show_ode(ftz,[a b],ya,50,'euler','b') %% these calculations show that, as the number of steps of the given interval %% increase, Euler's method converges to the solution, slowly but surely. %% I tried to squeeze a quick discussion of SHOOTING into the remaining 10mins % of class, but did not finish it, so, will pick it up next time.