[Content_Types].xml
_rels/.rels
matlab/document.xml
matlab/output.xml
metadata/coreProperties.xml
metadata/mwcoreProperties.xml
metadata/mwcorePropertiesExtension.xml
metadata/mwcorePropertiesReleaseInfo.xml
Stability Analysis In this livescript, you will learn how To apply stability analysis to determine the maximal time step that can be used to solve an initial value problem. Consider the linear initial value problem \frac{dx}{dt}=-100x with the initial condition x(0)=1 . Since the equation is linear, it’s quite easy to solve, especially if we use \texttt{dsolve} . syms x(t) t
eqn = diff(x,t) == -100*x;
cond = x(0) ==1;
sol = dsolve(eqn,cond) However, if we were to use a numerical time stepping method, such as the Explicit Euler method, the solution with a time step of \Delta t=0.02 is given by f = @(t,x) -100*x;
tspan = [0,0.1];
dt = 0.025;
x0 = 1;
t = tspan(1):dt:tspan(2);
x = zeros(length(t),1);
x(1) = x0;
for l = 1:length(t)-1
x(l+1,:) = x(l,:)+dt*f(t(l),x(l,:))’;
end Plotting against the analytical solution, we have x_exact = @(t) exp(-100*t);
fplot(x_exact,[0,0.1])
hold on
plot(tspan(1):dt:tspan(2),x)
hold off we can see that the solution obtained from the explicit Euler method diverges rapidly, which we call blow up. (a) Run the code with a reduced step size of \Delta t=0.02 and \Delta = 0.01 . So the question is: Why does the solution improve after we reduce the time step to \Delta t=0.02 . To see why this happens, let’s look at what happens after a single time step when \Delta t=0.025 . The exact solution is given by x_exact(0.025) which means that for the second step, we should be solving the initial value problem \frac{dx}{dt}=-100x with the initial condition x(0.025)=0.3679 . But for the explicit Euler method (this code is just the previous code with a time step of 0.025) x = explicit_euler_example(); What we can see is that the explicit Euler method deviates from the real solution by some amount, which we’ll call \epsilon . eps = x(2)-x_exact(0.025) So after one step with the explicit Euler method, we’re now solving a different initial value problem for the approximate position x_{a} , given by \frac{dx}{dt}=-100x with the initial condition x_{a}(0.025)=0.3679+\epsilon . Defining the error as e = x_{a}-x Subtracting the two differential equations, we find that the error evolves as \frac{de}{dt}=-100e with the initial condition e(0.025)=\epsilon . Using the results from lectures for the explicit Euler method, these errors will be damped out over time if the time step satisfies the inequality |1+\lambda h|\leq1 In our case, we have \lambda=-100 , which requires quite a small time step to remain stable. (d) Assuming that \lambda is real valued, show that the maximum time step is given by h\leq\frac{2}{|\lambda|} Hence, the maximum time step will be given by max_h = min(2./abs(-100)) which is what we found before. This shows us that the reason we look at the model problem \frac{dx}{dt}=\lambda x is because using a numerical method will always introduce some error, which peturbs the original equation to a different equation, hopefully still close by. But if the problem is poorly posed, then perturbing the problem can lead to different solutions. So, we therefore need to choose a time-step such that our numerical method will dampen out the errors, which is the aim of stability analysis. Note that stability and accuracy are different concepts, stability relates to how errors evolve over time, whereas accuracy is how well the approximation represents the exact solution. Nonlinear Problems Suppose now we complicate matters further by making the right hand side of the ODE nonlinear, such that \frac{dx}{dt}=f(t,x) Then applying a similar analysis as before, we can show that the error evolves as \frac{de}{dt}=f(t,x_{a})-f(t,x) This might look horrible to deal with, but if we take the Taylor series expansion of f(t,x_{a}) about the exact solution x , we find that \frac{de}{dt}=f(t,x)+\Big(\frac{\partial f}{\partial x}\Big)_{x}e+\mathcal{O}(e^{2})-f(t,x) Hence, to linear order, we have \frac{de}{dt}=\Big(\frac{\partial f}{\partial x}\Big)_{x}e Note that the derivative is time-dependent, so to determine \lambda=\partial f/\partial x , you’ll either have to determine upper and lower bounds for the derivative (which is quite tedious), or look at the physics of the problem to approximate \lambda at the point in time where it’s most likely to be an issue. function x = explicit_euler_example()
f = @(t,x) -100*x;
tspan = [0,1];
dt = 0.025;
x0 = 1;
t = tspan(1):dt:tspan(2);
x = zeros(length(t),1);
x(1) = x0;
for l = 1:length(t)-1
x(l+1) = x(l)+dt*f(t(l),x(l))’;
end
end
manual code ready 0.4 true false 0 8 8 false false 1 9 9 false false 2 10 10 false true 3 11 11 true false 4 14 14 false false 5 15 15 false false 6 16 16 false false 7 17 17 false false 8 19 19 false false 9 20 20 false false 10 21 21 false false 11 23 25 false false 12 27 27 false false 13 28 28 false false 14 29 29 false false 15 30 30 false true 16 31 31 true true 17 38 38 true false 18 44 44 false true 19 46 46 true true 20 63 63
2020-10-01T05:57:51Z 2020-10-11T02:15:22Z
application/vnd.mathworks.matlab.code MATLAB Code R2020b
9.9.0.1444674 13ec83af-3c87-4ee7-a421-6ad6cb0b8502
9.9.0.1467703
R2020b
Aug 26 2020
2314982500