IT1005

1. IT1005 Lab session on week 11 (7th meeting) 1 or 2 more weeks to go…

2. Lab 7 – Quick Check • Have you all received my reply for lab 7? • Of course not. I have not finished grading your submissions yet >.< • I have important research paper deadline >.< • http://www.cs.mu.oz.au/cp2008/ • Abstract due: 4 April 2008, Paper due: 8 April 2008. • My future reply should contains: • Remarks on your M files (look for SH7 tags again) • For CSTR_input.m, CSTR_matrix.m, and CSTR_soln.m, etc • Remarks on your Microsoft Word file • For the other stuffs • Your marks is stored in the file “Marks.txt” inside the returned zip file! • I do not think I will be strict mode , the answers are standard for L7! • Note that my marking scheme is slightly different from the standard one,as I put emphasize on coding style (indentation, white spaces), efficiency,things like proper plotting, etc…

3. L7.Q1 – 4 tanks CSTRs Q1. A. Transform the non standard set of Linear Equations into standard format: -(Q+k*V) *CA1 + 0 *CA2 + 0 *CA3 + 0 *CA4 = -Q*CA0 Q *CA1 -(Q+k*V) *CA2 + 0 *CA3 + 0 *CA4 = 0 0 *CA1 + Q *CA2 -(Q+k*V) *CA3 + 0 *CA4 = 0 0 *CA1 + 0 *CA2 + Q *CA3 -(Q+k*V) *CA4 = 0 Q1. B. Convert the standard format into matrix format. Another straightforward task: [ -(Q+k*V) 000 ] *[CA1] = [-Q*CA0 ] [ Q -(Q+k*V) 00 ] [CA2] [0 ] [ 0 Q -(Q+k*V)0 ] [CA3] [0 ] [ 0 0 Q-(Q+k*V) ] [CA4] [0 ] Q1. C. You just need to do: V = 1; Q = 0.01; CA0 = 10; k = 0.01; A = [-(Q+k*V) 0 0 0; Q -(Q+k*V) 0 0; 0 Q -(Q+k*V) 0; 0 0 Q -(Q+k*V)]; b = [-Q*CA0; 0; 0; 0]; x = A\b % you should get x = [5; 2.5; 1.25; 0.625], x = inv(A) * b is NOT encouraged! % Using fsolve for this is also not encouraged!

4. L7.Q2 – n tanks CSTRs Q2. A. CSTR_input.m Answer is very generic, similar to read_input.m in L6.Q2 (Car Simulation) Q2. B. CSTR_matrix.m % This is my geek version, do NOT USE THIS VERSION (too confusing for novice)! function [A b] = CSTR_matrix(Q, CA0, V, k, N) A = diag(repmat(-(Q+k*V),1,N)) + diag(repmat(Q,1,N-1),-1); % This is a bit crazy :$ b = zeros(N,1); b(1) = -Q*CA0; Q2. C. CSTR_soln.m clear; clc; clf; % New trick, but important ! Clear everything before starting our program! [Q CA0 V k N] = CSTR_input(); [A b] = CSTR_matrix(Q, CA0, V, k, N); plot(A\b,'o'); % I prefer not to connect the plot with line, but it is ok if you do so. title('CA of each tank'); xlabel('Tank no k'); ylabel('CA_k'); % Good for CSTR_plot.m axis([0.5 N+0.5 0 CA0]); % Fix y axis so that it is consistent across 4 plots (same CA0!)

5. L7.Q2 – Good Plot Remember: Plot A against B means thatA is the Y axis, B is the X axis! Should just stop here (n = 4) for case 1

6. L7.Q3 and Q4 • Tips for guessing logically: • Plot the functions with some range, see the region of zero intercepts… • Guess from the easiest equation!e.g. x + y + z = 5 • log(z)… hm… z should be > 0 Q3. Type these at command window: syms x y; % no need to say syms f1 f2, the next two lines will create f1 and f2 anyway f1 = x^2 * y^2 - 2*x - 1; f2 = x^2 - y^2 - 1; [a b] = solve(f1,f2); % or >> s = solve(f1,f2); a = s.x; b = s.y; eval(a), eval(b) % convert the symbolic values to numeric values, these are the roots Q4. Create this function function F = lab07d(x) F(1) = sin(x(1)) + x(2)^2 + log(x(3)) - 7; F(2) = 3*x(1) + 2^(x(2)) + 1 - x(3)^3; F(3) = x(1) + x(2) + x(3) - 5; At command window (wild guesses will likely give you many imaginary numbers): fsolve(@lab07d, [1 1 1])  x=0.5991, y=2.3959, z=2.0050 fsolve(@lab07d, [5 -1 1])  x = 5.1004, y = -2.6442, z= 2.5438

7. Application 5: IVP (revisited) • Equation • A statement showing the equality of two expressions usually separated by left and right signs and joined by an equals sign. • Differential Equation • A description of how something continuously changes over time. • Ordinary Differential Equation • A relation that contains functions of only one independent variable, and one or more of its derivatives with respect to that variable. • Initial Value Problem • An ODE together with specified value, called the initial condition, of the unknown function at a given point in the domain of the solution.

8. Application 5: IVP (revisited) • We have seen several examples of IVP throughout IT1005 • Spidey Fall example (Lecture 2, Lecture ODE 1) • How velocity v change over time dt? (dv/dt) • Depends on gravity minus drag! • How displacement s change over time dt? (ds/dt) • Depends on the velocity at that time. • The Initial Values for v and s at time 0  v(0) = 0; s(0) = 0; • Car accelerates up an incline example (Lab 6, Q2) • How velocity v change over time dt? (dv/dt) • Depends on engine force on some kg car minus friction and gravity factor! • How displacement s change over time dt? (ds/dt) • Depends on the velocity at that time. • The Initial Values for v and s at time 0  v(0) = 0; s(0) = 0; • And two more for Term Assignment (Q2 and Q3) 

9. Application 5: IVP (revisited) • Solving IVP (either 1 ODE or set of coupled ODEs): • Hard way/Traditional way/Euler method: • Time is chopped into delta_time, then starting from the initial values for each variable, simulate its changes over time using the specified differential equation! • What Colin has shown in Lecture 2 for Spidey Fall is a kind of Euler method. • What you have written for Lab 6 Question 2 (Car Simulation) is also Euler method. • Matlab IVP solvers: mostly numerical solutions for ODE. • Create a derivative function to tell Matlab about how a variable change over time! function dydt = bla(t,y) % always have at least these two arguments % explain to Matlab how to derive dy/dt! Can be for coupled ODEs! dydt = dydt'; % always return a column vector! • Call one of the ODE solver with certain time span and initial values [t, y] = ode45(@bla, [tStart tEnd], IVs); % IVs is a column vector for coupled ODEs! plot(t,y(:,1)); % we can immediately plot the results (also in column vector!)

10. Term Assignment – Admin • This is 30% of your final IT1005 grade... Be very serious with it. • No plagiarism please! • Even though you can ‘cross check’ your answers with your friends (we cannotprevent that), you must give a very strong ‘individual flavor’ in your answers! • The grader will likely grade number per number, so he will be very curiousif he see similar answers across many students. Do not compromise your 30%! • Who will grade our term assignment? • I may not be the one doing the grading! Perhaps all the full time staff… dunno yet. • Submit your zip file (containing all files that you use to answer the questions)to IVLE “Term Assignment” folder! NOT to my Gmail! • Your zip file name should be: yyy-uxxxxxx.zip, NOT according to my style! • Strict deadline, Saturday, 5 April 08, 5pmThat IVLE folder will auto close by Saturday, 5 April 08, 5.01pmBe careful with NETWORK CONGESTION around these final minutes…To avoid that problem, submit early, e.g. Friday, 4 April 08, night.

11. Term Assignment – Q1 • Question 1: Trapezium rule for finding integration • A. Naïve one. Explain your results! • B. More accurate one. Explain your results! • References: • help quad • http://en.wikipedia.org/wiki/Numerical_integration • http://en.wikipedia.org/wiki/Trapezium_rule • Revision(s) to the question: • Symbol ‘a’ changed to ‘c’ inside function f(t)! • In Q1.B, the rows in column ‘c’ are [0.001 0.5 10.0 100.0] not [0.01 0.5 1.0 10.0]! • In Q1.B, the range of k is changed from k = 2:n-1 to k = 1:n-1,but it is ok if you stick with the old one!

12. Term Assignment - Q2 • Question 2: Zebra Population versus Lion Population • A. IVP, coupled, non-linear ODEs • B. Explain what you see in the graph of part A above. • C. Steady state issue. • D. IVP again, but change the IVs according to part C above. Comment! • E. IVP with different IVs, and different plotting method. Comment! • References: • Google the term ‘predator prey’ as mentioned in the question. • help odeXX (depends on the chosen solver) • http://en.wikipedia.org/wiki/Steady_state • Revision(s) to the question: • No change so far…

13. Term Assignment – Q3 • Question 3: Similar to Q2, Predator-Prey: n = 4 species • A. IVP again, 4 coupled, non-linear ODEs. Dr Saif said that we must use ode15s! (See ODE 3 & 4 lecture note) • B. IVP, same IVs, 1.000 years, 3D plot x1-x2-x3 (x4 is not compulsory),and explain. • C. Explain plot in B as best as you can. • References: • http://en.wikipedia.org/wiki/Lotka-Volterra_equation (mentioned in the question). • Google ‘Matlab 3D plot’ • help odeXX (depends on the chosen solver) • Revision(s) to the question: • The ODE equations are updated! Read the newest one! • The coefficient r(3) is changed from 1.53 to 1.43!

14. Free and Easy Time • Now, you are free to explore Matlab to: • Do your Term Assignment (all q1, q2, and q3) • You should NOT use me as an oracle, e.g. • I cannot find the bug in my program, can you help me? • Are my 2-D/3-D plots correct? • Are my …. bla bla … correct?