MATLAB Program - Assignment Example

EG-228 EG-264 I confirm that I have not received help from, or given help to, anyone else in constructing the solution to this assignment. here and signature above) 530542 1. The force, F(x), in the spring when it is compressed by an amount x m by the mass is given to be F(x) = k1xa + k2xb Where k1, k2, a and b are constants. Produce a plot, showing the variation of F with x, for 0 ≤ x ≤ 6h. To answer this, we must first compute the values of k1, k2, a, b and h. k1 = 45000 (1 + u) = 45000 (1 + 5) = 270000 k2 = 8100 (1 + v) = 45000 (1 + 3) = 32400 a = 1 + 0.01(z) = 1 + 0.01(2) = 1.02 b = 1.5 + 0.01(u) = 1.5 + 0.01(5) = 1.55 h = 0.43 (1 + 0.1*x) = 0.43 (1 + 0.1*5) = 0.645 Thus, the equation we want to graph is F(x) = 270000x1.02 + 32400x1.55, for values of x between 0 and 6*0.645. The following MATLAB program produces the plot: function plot1() X = linspace(0,6*0.645,1000); Y = f(X); plot(X,Y); end function [f_value]=f(x) f_value = 270000*x.^(1.07) + 32400*x.^(1.55); end The book by Steven Karris, although intended for students in Digital Signal Processing course, provides an excellent tutorial for MATLAB beginners in its appendices. Karris (A-14) explains that the format for linspace command is linspace(first_value, last_value, number_of_values). We use this command to define the upper and lower limits of the graph, as well as the number of values sufficient to produce a smooth graph. The third line calls the subfunction f_value where the function to be graphed is defined. The following figure shows the plot produced: 2. Conservation of energy can be used to show that, at equilibrium, where d is the compression of the spring. Use the bisection method to determine an approximation to the value of d, with an absolute error of less than 1 x 10-6. First, we solve for the necessary variables: m = 960 (1.2 – 0.1w) = 960 (1.2 – 0.1*0) = 1152 g = 9.8 (1 + 0.01y) = 9.8 (1 + 0.01*4) = 10.192 For this problem, we need not write a new MATLAB program. We simply reuse the bisection program we’ve used before. However, we cannot directly plug in the equation as it needs to be expressed in the form F(x) = 0. Thus the equation becomes . function bisection() clc n = 21; % specify number of bisection steps % define the initial interval a = 0; b = 2; fprintf(\n initial interval [%g, %g] \n total bisection steps %d\n, a,b,n); x_left = a; x_right = b; f_left = f(x_left); f_right = f(x_right); if f_left*f_right > 0 error(ERROR: no root in the specified interval); end %the bisection method for i=1:n if f_left == 0 %exact root is reached by the left bound of the interval fprintf(\n stage %g root %g with zero absolute error \n,i,x_left); return; end if f_right==0 %exact root is reached by the right bound of the interval fprintf(\n stage %g root %g with zero absolute error \n,i,x_right); return end %the bisection process x_mid = (x_left+x_right)/2.0; f_mid = f(x_mid); if f_left*f_mid 0 % there is a root in [x_mid,x_right] x_left = x_mid; f_left = f_mid; end %compute the approximate root for the current bisection step root = (x_left+x_right)/2.0; %compute the absolute error for the current bisection step abs_err=(x_right-x_left)/2.0; fprintf(\n stage %g root %g absolute error < %g \n,i,root,abs_err); end %check satisfaction of equation at end of process residual = f(root); fprintf(\n final residual = %g \n,residual); end %Subfunction defines the equation f(x) = 0 function f_value = f(x) f_value = (270000/2.02)*x.^(2.02) + (32400/2.55)*x.^(2.55) - 1152*10.192*x - 1152*10.192*0.645; end As can be seen, not much has changed in our program except in the subfunction that defines the equation and the values of n, a, and b. The number of bisection steps that gives an error less than 1 x 10-6 is 21 when the interval is 0 to 2. The actual output is as follows: initial interval [0, 2] total bisection steps 21 stage 1 root 0.5 absolute error < 0.5 stage 2 root 0.25 absolute error < 0.25 stage 3 root 0.375 absolute error < 0.125 stage 4 root 0.3125 absolute error < 0.0625 stage 5 root 0.28125 absolute error < 0.03125 stage 6 root 0.296875 absolute error < 0.015625 stage 7 root 0.289063 absolute error < 0.0078125 stage 8 root 0.285156 absolute error < 0.00390625 stage 9 root 0.283203 absolute error < 0.00195313 stage 10 root 0.282227 absolute error < 0.000976563 stage 11 root 0.282715 absolute error < 0.000488281 stage 12 root 0.282471 absolute error < 0.000244141 stage 13 root 0.282349 absolute error < 0.00012207 stage 14 root 0.282288 absolute error < 6.10352e-005 stage 15 root 0.282257 absolute error < 3.05176e-005 stage 16 root 0.282242 absolute error < 1.52588e-005 stage 17 root 0.282249 absolute error < 7.62939e-006 stage 18 root 0.282246 absolute error < 3.8147e-006 stage 19 root 0.282248 absolute error < 1.90735e-006 stage 20 root 0.282248 absolute error < 9.53674e-007 stage 21 root 0.282249 absolute error < 4.76837e-007 final residual = 0.00633197 3. The energy stored in the spring, at equilibrium, may be expressed in the form Use the composite trapezoidal rule to obtain an approximation to the value of E, with a relative error of less than 0.001%. Once again, we can use the existing composite trapezoidal program for this problem. Unlike problem 2 however, a considerable amount of codes has to be added as the program should be able to automatically determine the number of intervals that will give a relative error less than 0.001 percent. This feature was not available in the original program. Modified program: function trapezoidal(integ) %INTTRAPEZOIDAL integration using the composite trapezoidal rule clc a = 0; b = 0.282249; if (nargin ~= 1) integ = quad(@f, a, b); end n = 1; rel_err = 100; err = 0.001; while (err < rel_err) integral = 0.0; h = (b-a)/n; for m = 1:n x_left = a+(m-1)*h; x_right = a+m*h; f_left = f(x_left); f_right = f(x_right); integral = integral+h/2*(f_left+f_right); end rel_err = (abs(integral – integ)/integ)*100; n = n+1; end fprintf(‘\n Trapezoidal approximation to integral using %g subintervals is %g \n’,n-1,integral); fprintf(‘\n Relative error is %g percent.\n’, rel_err); end function f_value = f(x) f_value = 270000*x.^(1.07) + 32400*x.^(1.55); end Because the problem requires the calculation of percent relative error, we need to know the exact answer to the integral. Percent relative error is defined as The program was designed such that it accepts the exact answer from the user’s input (the user should call it “integ”) or calculates the answer using the quad command when there is no input coming from the user. The quad command performs numerical integration using Simpson’s Rule (8-59). Next, we declare the initial value of relative error as 100 percent. Since we have no clue as to what number of intervals give a relative error of less than 0.001 percent, we initialize the value of n as 1. The while loop keeps the program iterating starting from n = 1 and incrementing n by 1 with each iteration until the percent relative error is less than the prescribed relative error (0.001 percent). The output is as follows if we input 10015.3 as the exact answer: >> trapezoidal(10015.3) Trapezoidal approximation to integral using 86 subintervals is 10015.4 Relative error is 0.000978298 percent. Meanwhile, the following is the output when the user does not give the exact answer: Trapezoidal approximation to integral using 94 subintervals is 10015.4 Relative error is 0.000997879 percent. As can be seen, both give the same answer but different values of calculated subintervals and percent relative error. Works Cited Karris, Steven T. Signals and Systems with MATLAB Computing and Simulink Modeling, 3rd ed. USA: Orchard Publications, 2007. Read More