channelflow.org

====== Differences ====== This shows you the differences between two versions of the page.

--- gibson:teaching:spring-2016:math445:lecture:timestepping [2016/04/14 09:53]
gibson
+++ gibson:teaching:spring-2016:math445:lecture:timestepping [2016/04/14 18:23] (current)
gibson [Problem 4]
@@ Line 79: / Line 79: @@
 </file>
- {{ :gibson:teaching:spring-2016:math445:lecture:cylinderpath2.png?direct&400 }}
+ {{ :gibson:teaching:spring-2016:math445:lecture:cylinderpath1.png?direct&400 }}
+Note that the trajectory computed here is not very accurate. The particle shouldexit the box at the same $y$ value it had when it entered. The problem is we chose quite a large time step $\Delta t = 0.4$, and forward-Euler is only 1st-order accurate (error scales as $\Delta t$). In the next problem, we'll reduce the time step to $\Delta t = 0.01$ and get a more accurate solution --though still not as good as the 4th-order accurate ''ode45'' function.
-Note that the computed trajectory is not very accurate, since we chose quite a large time step $\Delta t = 0.4$, and forward-Euler is only 1st-order accurate (error scales as $\Delta t$).
 ----
 ====Problem 2====
-Write Matlab code that plots the //path// of the particle as a red curved line. To do this wee need to save the sequence of $\vec{x}$ values in a matrix, and then plot the rows of that matrix as a line.
+Write Matlab code that plots the //path// of the particle as a red curved line. To do this we need to save the sequence of $\vec{x}$ values in a matrix, and then plot the rows of that matrix as a line.
 <code matlab>
@@ Line 107: / Line 109: @@
 ylim([-2,2])
 </code>
+ {{ :gibson:teaching:spring-2016:math445:lecture:cylinderpath3.png?400 |}}
+----
+====Problem 3====
+Write Matlab code that plots the path of the particle as a red curved line, using Matlab's ''ode45'' function to do the time-integration.
+<code matlab>
+% use Matlab's ode45 function to do the time integration of
+% dx/dt = v(t, x)
+T = 10;           % integrate from t=0 to T=10
+x0 = [-2.8; 0.6]; % initial position of particle
+[t, x] = ode45(@v, [0 T], x0);  % computes x(t) at given values of t
+plot(x(:,1), x(:,2), 'r-');
+axis equal
+axis tight
+xlim([-3,3])
+ylim([-2,2])
+</code>
+Matlab's ''ode45'' function requires an ODE of the form $dx/dt = v(t,x)$, so we have to add a ''t'' argument to our ''v'' function
+<file matlab v.m>
+function dxdt = v(t, x);
+  % compute 2D inviscid cylinder velocity as function of x
+  % input: vector x has x(1) = x coord, x(2) = y coord
+  V0 = 1;  % scale of velocity
+  a  = 1;  % radius of circle
+  % compute polar coordinates
+  r = sqrt(x(1)^2 + x(2)^2);
+  theta = atan2(x(2), x(1));
+  dxdt = [V0*(1 - (a/r)^2 * cos(2*theta)) ;
+         -V0*(a/r)^2 * sin(2*theta)];
+end
+</file>
+This produces a plot very like the one for Problem 3.
+----
+==== Problem 4 ====
+Draw a number of particle paths starting with a number of $y$ values and $x=-2.8$.
+<code matlab>
+% use Matlab's ode45 function to do the time integration of
+% dx/dt = v(t, x)
+T = 10;                           % integrate from t=0 to T=10
+for y = -2:0.4:2                  % loop over different y values
+  x0 = [-2.8; y];                 % set initial position of particle
+  [t, x] = ode45(@v, [0 T], x0);  % compute x(t) over range 0 <= t <= T
+  plot(x(:,1), x(:,2), 'r-');     % plot the path
+end
+</code>
+{{ :gibson:teaching:spring-2016:math445:lecture:cylinderpath2.png?direct&400 |}}

channelflow.org

User Tools

Site Tools

Page Tools