Unformatted text preview: 1/14/09 The most basic operator for making decisions is the if statement. We will use it to implement the bisection method of Section 5.4. In this method we start with a function that is positive at x=xl and negative at x=xu (or vice versa) and we evaluate it in the middle. Then depending on the sign, we repeat the same procedure with either the left or right interval. We will start with the assumption that the function is positive at xl, and develop an Mfile for performing one step of the bisection for the cosine function. function [xl,xu]= cosbisect(xl,xu) %performs one step of bisection towards finding the root of cos(x) in the %interval (xl,xu) xm=(xl+xu)/2; if (cos(xm) < 0) xu=xm; else xl=xm; end Then we can invoke it repeatedly to narrow down the interval >> >> xl=0;xu=2; >> [xl,xu]= cosbisect(xl,xu) xl = 1.5000 xu = 2 >> [xl,xu]= cosbisect(xl,xu) xl = 1.5000 xu = 1.7500 >> [xl,xu]= cosbisect(xl,xu) xl = 1.5000 xu = 1.6250 >> [xl,xu]= cosbisect(xl,xu) xl = 1.5625 xu = 1.6250 >> [xl,xu]= cosbisect(xl,xu) xl = 1.5625 xu = 1.5938 Note that if we do not call the function with [xl,xu] on the right side, it will not update them >> xl=0;xu=2; >> cosbisect(xl,xu) ans = 1 >> xl,xu xl = 0 xu = 2 Instead of the if statement, we can also use the switch statement. The latter is useful when there are more than two possibilities based on a single test quantity. In the previous example, we excluded the possibility that the cosine of the middle point will be zero. With the switch structure we can cater to that possibility and tell the user that the root has been found. function [xl,xu]= cosbisect_case(xl,xu) %performs one step of bisection towards finding the root of cos(x) in the %interval (xl,xu) xm=(xl+xu)/2; switch sign(cos(xm)) case 1 xu=xm; case 1 xl=xm; case 0 disp('root found') disp(xm) end And we can exercise it as before: >> [xl,xu]= cosbisect_case(xl,xu) xl = 1 xu = 2 >> [xl,xu]= cosbisect_case(xl,xu) xl = 1.5000 xu = 2 >> [xl,xu]= cosbisect_case(xl,xu) xl = 1.5000 xu = 1.7500 >> [xl,xu]= cosbisect_case(xl,xu) xl = 1.5000 xu = 1.6250 If we want to carry this operation for a specified number of times we can put it in a loop. The simplest loop in Matlab is the for loop where the number of repetition is specified in the form for index start:step:finish end. statements So we can change the bisection function to execute for a given number of times and store the history of the middle point. function [xl,xu,xm,ym]= cosbisect_loop(xl,xu,n) %performs n steps of bisection towards finding the root of cos(x) in the %interval (xl,xu), updates the interval, and stores the values of the %midpoint and the function for i=1:1:n xm(i)=(xl+xu)/2; ym(i)=cos(xm(i)); switch sign(ym(i)) case 1 xu=xm(i); case 1 xl=xm(i); case 0 disp('root found') disp(xm) end end Then we exercise it to get [xl,xu,xm,ym]=cosbisect_loop(0,2,6) xl = 1.5625 xu = 1.5938 xm = 1.0000 1.5000 1.7500 1.6250 1.5625 1.5938 ym = 0.5403 0.0707 0.1782 0.0542 0.0083 0.0230 >> [xl,xu,xm,ym]=cosbisect_loop(0,2,10) xl = 1.5703 xu = 1.5723 xm = 1.0000 1.5000 1.7500 1.6250 1.5625 1.5938 1.5781 1.5703 1.5742 1.5723 ym = 0.5403 0.0707 0.1782 0.0542 0.0083 0.0230 0.0073 0.0005 0.0034 0.0015 We can also call the function with a reduced set of arguments [xl,xu]=cosbisect_loop(0,2,10) xl = 1.5703 xu = 1.5723 >> Instead of a loop with a given number of repetitions, it may make more sense to exit when a condition is satisfied. This can be done with the while structure. For finding a root, it is natural to condition the loop on the accuracy of the solution. In the case of the bisection method, it is the size of the interval function [xm,ym]= cosbisect_while(xl,xu,accuracy) %performs bisection towards finding the root of cos(x) in the %interval (xl,xu), until the interval is smaller than accuracy and outputs the value of the %midpoint and the function while (xuxl)>accuracy xm=(xl+xu)/2; ym=cos(xm); if sign(ym)<0 xu=xm; else xl=xm; end end Exercising with several values of the accuracy: >> [xm,ym]=cosbisect_while(0,2,0.1) xm = 1.5625 ym = 0.0083 >> [xm,ym]=cosbisect_while(0,2,0.01) xm = 1.5703 ym = 4.8383e004 >> [xm,ym]=cosbisect_while(0,2,0.001) xm = 1.5713 ym = 4.9274e004 >> [xm,ym]=cosbisect_while(0,2,0.0001) xm = 1.5707 ym = 5.6581e005 ...
View
Full Document
 Spring '09
 RAPHAELHAFTKA
 xl,xu

Click to edit the document details