Question: Using your preferred code editor (e.g. VSCode), in a Python script called CR3.py , write code to answer the following problem. Problem Consider the
Using your preferred code editor (e.g. VSCode), in a Python script called "CR3.py", write code to answer the following problem.
Problem
Consider the following ODE, describing the displacement u(t) of an oscillator (e.g. a mass attached to a spring) at time t0 away from its resting position:
d^2u/dt^2=((0)^2)*u,
where 0 is the angular frequency of the oscillator, given as a fixed parameter. The initial conditions are given as
- u(t=0)=u0, the initial displacement,
- du/dt(t=0) = v0, the initial velocity.
We seek to solve this equation numerically, using a finite difference method. We discretise time with a step size t, and denote by Un the approximation to u(t) computed at time step n, that is Unu(t=n*t)
We approximate the second derivative in the ODE with DC^2(t), to obtain the difference equation:
(U(n+1)2*Un+U(n1))/(t^2)=(0)^2)*Un.
Rearranging this equation allows us to compute Un+1 using values computed at the two previous time steps, Un and Un1. We can initialise U0 and U1 by discretising the initial conditions:
- U0=u0,
- (U1U0)/t=v0.
Then, we can compute U2, followed by U3, etc., for as many time steps as we desire.
Note that this method provides a "valid" solution as long as 0*t<2. This is called the stability condition.
Your task
Write a function "oscillator(w0, u0, v0, nmax, dt)" which takes 5 input arguments:
- a positive number "w0", representing the parameter 0,
- two numbers "u0" and "v0", representing the initial conditions u0 and v0 respectively,
- a positive integer "nmax" greater than or equal to 2, representing the total number of time steps,
- a positive number "dt", representing the step size t,
and returns a Numpy vector "U" with a total of "nmax" elements, where the iith element is the value of Ui, the approximated solution for the oscillator displacement at the iith time step, computed using the method described above.
Additionally, before computing the result, the function should check the value of t and display a message if necessary:
- if 0*t>2, display a message warning the user that the stability condition has been violated. Your message should include the maximum value of "dt" allowed for the chosen value of "w0", to help your user choose a better step size next time.
- if 0*t=2, display a message warning the user that the chosen step size is at the stability limit.
In either case, you should still compute and return the solution.
Testing
After the function definition, write a few tests to check that your function is working. Your first test should use the following values:
0=5,u0=0.2,v0=5,nmax=500,t=0.03.
You can compare the computed solution to the exact solution, which is given by:
u(t)=u0*cos(0*t)+(v0/0)*sin(0*t).
You could plot the exact and computed solutions on the same graph over time -- for valid step sizes (such that 0*t<2), the two curves should be close, but they will never be exactly overlapping; the peak values will often be where you see the largest error. For values of t and 0 close to the stability condition (but still valid), you should see that the computed solution oscillates slightly faster than the exact solution. As we have seen in many other situations, decreasing the step size t should generally improve accuracy.
Try different values for the initial conditions u0 and v0, and check that the first two values of your computed solution are set correctly. Try different values for the frequency 0, increasing it should produce a faster-oscillating solution.
Make sure you check that the warning messages appear for appropriate values. Plot the computed solution for both problematic cases -- this is what you should see:
- when (and only when) 0*t>2, the amplitude of the oscillations in the computed solution should grow exponentially over time, and you should see the first warning message.
- when (and only when) 0*t=2, the amplitude of the oscillations in the computed solution should grow linearly over time, and you should see the second warning message.
- when (and only when) 0*t<2, the computed solution should oscillate without growing linearly or exponentially over time, and you should not see any message.
Step by Step Solution
There are 3 Steps involved in it
Get step-by-step solutions from verified subject matter experts
