gibson:teaching:spring-2016:math445:lecture:timestepping
====== Differences ====== This shows you the differences between two versions of the page.
| Both sides previous revision Previous revision Next revision | Previous revision | ||
|
gibson:teaching:spring-2016:math445:lecture:timestepping [2016/04/14 09:46] gibson |
gibson:teaching:spring-2016:math445:lecture:timestepping [2016/04/14 18:23] (current) gibson [Problem 4] |
||
|---|---|---|---|
| Line 11: | Line 11: | ||
| The above equation is an [[https://en.wikipedia.org/wiki/Ordinary_differential_equation| ordinary differential equation]]. MATH 527 //Differential Equations// is all about finding solutions to such equations in closed form. Here we will use **numerical solution methods** to find approximate numerical solutions. | The above equation is an [[https://en.wikipedia.org/wiki/Ordinary_differential_equation| ordinary differential equation]]. MATH 527 //Differential Equations// is all about finding solutions to such equations in closed form. Here we will use **numerical solution methods** to find approximate numerical solutions. | ||
| - | To make things concrete, let's solve the following problem. Given the 2D inviscid flow around a cylinder from the [[gibson:teaching:spring-2016:math445:lecture:cylinderflow|last lecture]], what would be the path of a small particle starting at position $(x,y) = (-3.8, 0.6)$, shown below as a red circle? | + | To make things concrete, let's solve the following problem. Given the 2D inviscid flow around a cylinder from the [[gibson:teaching:spring-2016:math445:lecture:cylinderflow|last lecture]], what would be the path of a small particle starting at position $(x,y) = (-2.8, 0.6)$, shown below as a red circle? |
| - | {{ :gibson:teaching:spring-2016:math445:lecture:cylinderpath1.png?direct&400 }} | + | {{ :gibson:teaching:spring-2016:math445:lecture:cylinderpath0.png?direct&400 }} |
| [[gibson:teaching:spring-2016:math445:lecture:cylinderflow_code|Here's the code for plotting the 2D inviscid cylinder flow]]. The line ''x = [-3.8, 0.6]; plot(x(1), x(2), 'ro');'' plots the dot. The problem now is to compute and plot the particle as it's advected along by the fluid velocity field $v(x)$, where | [[gibson:teaching:spring-2016:math445:lecture:cylinderflow_code|Here's the code for plotting the 2D inviscid cylinder flow]]. The line ''x = [-3.8, 0.6]; plot(x(1), x(2), 'ro');'' plots the dot. The problem now is to compute and plot the particle as it's advected along by the fluid velocity field $v(x)$, where | ||
| Line 52: | Line 52: | ||
| Nsteps = round(T/dt); | Nsteps = round(T/dt); | ||
| - | x = [-3.8; 0.6]; % initial position of particle | + | x = [-2.8; 0.6]; % initial position of particle |
| hold on; | hold on; | ||
| Line 79: | Line 79: | ||
| </file> | </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. | ||
| ---- | ---- | ||
| Line 85: | Line 88: | ||
| ====Problem 2==== | ====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> | <code matlab> | ||
| Line 106: | Line 109: | ||
| ylim([-2,2]) | ylim([-2,2]) | ||
| </code> | </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 |}} | ||
| + | |||
gibson/teaching/spring-2016/math445/lecture/timestepping.1460652360.txt.gz · Last modified: 2016/04/14 09:46 by gibson
Page Tools
Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Share Alike 3.0 Unported
