Homework #13

15.17

Matlab:

w1 = 2;

w2 = 2;

L = @(th,al) w1./sind(th)+w2./sind(180-al-th);

hold off;

L_thmin = [];

for al = [45:135],

th = [10:179-135];

plot(th,L(th,al));

hold on;

[lm,Lin] = min(L(th,al))

L_thmin = [L_thmin,min(L(th(Lin),al))];

end

hold off

plot((45:135),L_thmin);

16.32

For this problem the derivative was taken for the function then the value was places back in the original value to find the maximum.

Matlab:

e_0 = .85e-12;

q = 2e-5;

a = 0.9;

x = (0:0.001:15);

Q = q;

c1 = 1/(4*pi^e_0) *(q*Q);

F =@(x) c1*x./(x.^2+a.^2).^(3/2);

hold off;

plot (x,F(x));

hold on;

plot (a/sqrt(2), F(a/sqrt(2)), '*r');

16.25

The solution can be found at 1,3 based on using the max function in matlab.

Matlab:

s = (0:0.001:10);

T =@(s) 15* (s-s.^2)./((1-s).*(4*s.^2-3*s+4));

hold off

plot (s,T(s));

ess = rand(1,100000) *10;

[y,i] = max(T(ess));

hold on

plot (ess(i),y,'r*');

Homework # 12

14.5

Find the gradient vector and Hessian matrix for each of the following functions

(a) f(x,y) = 2xy^2 + 3e^(xy)

(b) f(x,y,z) = x^2 + y^2 + 2z^2

(c) f(x,y) = ln(x^2 + 2xy + 3y^2)

14.6

The returned minimum value was -0.5 for the function:

Here is the matlab code:

[X,Y] = meshgrid(0:.2:8, -4:.2:8);

f_xy = @(x,y) (x - 3).^2 + (y - 2).^2

df_xy_dx = @(x,y) 2*(x-3);

df_xy_dy = @(x,y) 2*(y-2);

x_0 = 0;

y_0 = 0;

g_h = @(h) f_xy(x_0+df_xy_dx(x_0,y_0)*h,y_0+df_xy_dy(x_0,y_0)*h)

figure(1);

h = -1:0.001:0;

hold off;

plot(h,g_h(h))

[gmin,hi] = min(g_h(h));

hold on 

plot(h(hi),gmin,'k*');

x_1 = x_0 + df_xy_dx(x_0,y_0)*h(hi)

y_1 = y_0 + df_xy_dy(x_0,y_0)*h(hi)

14.8

    This is the Matlab code we went over in class that shows the inital step along with two more steps to find the minimum for the given function

[X,Y] = meshgrid(-4:.2:8, -4:.2:8);

f_xy = @(x,y) -8*x+x.^2+12*y+4*y.^2-2.*x.*y

df_xy_dx = @(x,y) -8 + 2*x-x.*y;

df_xy_dy = @(x,y) 12+8*y-2*x;

figure(1)

hold off

surfc(X,Y,f_xy(X,Y))

hold on

x_0 = 0;

y_0 = 0;

g_h = @(h) f_xy(x_0+df_xy_dx(x_0,y_0)*h,y_0+df_xy_dy(x_0,y_0)*h)

figure(2);

h = -4:0.001:4;

hold off;

plot(h,g_h(h))

[gmin,hi] = min(g_h(h));

hold on 

plot(h(hi),gmin,'k*');

x_1 = x_0 + df_xy_dx(x_0,y_0)*h(hi)

y_1 = y_0 + df_xy_dy(x_0,y_0)*h(hi)

figure(1)

plot3(x_0,y_0,f_xy(x_0,y_0),'w*')

plot3(x_1,y_1,f_xy(x_1,y_1),'y*')

g_h = @(h) f_xy(x_1+df_xy_dx(x_1,y_1)*h,y_1+df_xy_dy(x_1,y_1)*h)

figure(3);

h = -4:0.001:4;

hold off;

plot(h,g_h(h))

[gmin,hi] = min(g_h(h));

hold on 

plot(h(hi),gmin,'k*');

x_2 = x_1 + df_xy_dx(x_1,y_1)*h(hi)

y_2 = y_1 + df_xy_dy(x_1,y_1)*h(hi)

figure(1)

plot3(x_1,y_1,f_xy(x_1,y_1),'y*')

plot3(x_2,y_2,f_xy(x_2,y_2),'W*')

g_h = @(h) f_xy(x_2+df_xy_dx(x_2,y_2)*h,y_2+df_xy_dy(x_2,y_2)*h)

figure(4);

h = -4:0.001:4;

hold off;

plot(h,g_h(h))

[gmin,hi] = min(g_h(h));

hold on 

plot(h(hi),gmin,'k*');

x_3 = x_2 + df_xy_dx(x_2,y_2)*h(hi)

y_3 = y_2 + df_xy_dy(x_2,y_2)*h(hi)

figure(1)

plot3(x_2,y_2,f_xy(x_2,y_2),'W*')

plot3(x_3,y_3,f_xy(x_3,y_3),'y*')

Homework #11

13.7

For this problem I first graphed he equation to see if there was a maximum for that range of x, and there was so I implemented the Golden-Section Search and reached the convergence of the true value of 0.8408 at x= -0.3472

Matlab:

%Graph that shows theres a maximum between -2:1

x = (-2:0.1:1);

f = @(x) -x.^4 - 2*x.^3 - 8*x.^2 - 5*x;

plot(x,f(x));

%plot(x,f_x(x));

xl = -2;

xu = 1;

 d = ((sqrt(5) - 1)/2)*(xu-xl);

 x1 = xl + d;

 x2 = xu - d;

 sol= [xl f(xl) x2 f(x2) x1 f(x1) xu f(xu) d]

 for iter = 1:13

 if f(x2) > f(x1)

    xu = x1;

    x1 = x2;

    d = ((sqrt(5) - 1)/2)*(xu-xl);

    x2 = xu - d;

 end

 if f(x1) > f(x2)

    xl = x2;

    x2 = x1;

    d = ((sqrt(5) - 1)/2)*(xu-xl);

    x1 = xl + d; 

 end

 sol= [xl f(xl) x2 f(x2) x1 f(x1) xu f(xu) ]

 end

 hold on;

 plot(x2,f(x2),'*r')

13.8

(a) for this problem I used the some code you provided us and altered it to work with this problem. 

Output:

GoldenRatio = -0.3472    0.8408

Parabolic = -0.3473    0.8408

Newton =  -0.3473    0.8408

Matlab:

theta0 = 50;

x = [-2:0.01:1];

y = @(x) -x.^4 - 2*x.^3 - 8*x.^2 - 5*x;

figure( 1);

plot(x,y(x));

hold on;

% Golden Ration Search

xl= -2;

xu = 1;

for ii = 0:26,

d = (xu-xl)*(sqrt(5)-1)/2;

x1 = xl+d;

x2 = xu-d;

if(y(x1)>y(x2))

xl = x2;

else

xu = x1;

end;

plot(xu,y(xu),'*r')

plot(xl,y(xl),'*r')

end

plot(xu,y(xu),'*g')

GoldenRatio = [xu y(xu)]

figure(2)

x = [0:0.01:20];

plot(x,y(x));

hold on;

x0 = -2;

x1 = -1;

x2 = 1;

% parabolic interpolation

for ii = 0:6,

x3 = y(x0)*(x1^2-x2^2) + y(x1)*(x2^2-x0^2)+y(x2)*(x0^2-x1^2);

x3 = x3/(2*y(x0)*(x1-x2)+2*y(x1)*(x2-x0)+2*y(x2)*(x0-x1));

if(abs(x3-x1)<1.0e-6) break; end

if(y(x3)>y(x1))

x0 = x1;

x1 = x3;

% x2 = x2;

else

x2 = x1;

x1 = x3;

% x0 = x0; 

end;

plot(x3,y(x3),'*r')

end

plot(x3,y(x3),'*g')

Parabolic = [x3 y(x3)]

% Newtons

figure(3)

y = @(x) -x.^4 - 2*x.^3 - 8*x.^2 - 5*x;

yp = @(x) -4*x.^3 -6*x.^2 -16*x -5;

ypp = @(x) -4*(3*x.^2 + 3*x +4);

x = [0:0.01:20];

plot(x,y(x),x,yp(x))

hold

x = 0;

for ii = 0:6,

x = x-yp(x)/ypp(x);

plot(x,y(x),'*r')

end

Newton = [x y(x)]

plot(x,y(x),'*g')

13.21

For this problem I implemented the golden section search to calculate the time and altitude of peak elevation

Output:

GoldenRatio = 3.8699(time)  192.8537(elevation)

Matlab:

g = 9.81;

z0 = 100;

v0 = 55;

m = 80;

c= 15;

t = (0:0.01:10);

f = @(t) z0 + (m/c)*(v0+((m*g)/c))*(1-exp(-((c/m)*t))) - ((m*g)/c)*t;

figure(1);

plot(t,f(t));

hold on;

% Golden Ration Search

xl= 0;

xu = 10;

for ii = 0:10,

d = (xu-xl)*(sqrt(5)-1)/2;

x1 = xl+d;

x2 = xu-d;

if(f(x1) > f(x2))

xl = x2;

else

xu = x1;

end;

plot(xu,f(xu),'*r')

plot(xl,f(xl),'*r')

end

plot(xu,f(xu),'*g')

GoldenRatio = [xu f(xu)]

13.22

Homework # 7

7.17

First I solved it graphically and I can approximate the root to 2.55

Matlab:

x = (-8:0.0001:6);

f_x = @(x) x.^3 + 3.5*x.^2 -40;

plot(x,f_x(x));

Then I solved it using the function fzero in Matlab with the result of 2.5676

Matlab:

x = fzero(inline('x.^3 + 3.5*x.^2 -40'),3)

 7.19

I used Wolfram Alpha to calculate the roots for the first part of the equation:

Then I plugged the roots in and got this equation:

8.9

The First thing I did was rearrange the equation to this:

Then I used the zero function in Matlab to solve the equation and I got: 2.9449

Matlab:

x = fzero(inline('(pi*h.^3)/3 - pi*(h.^2)*(1) + (0.5)'),3)

8.14

Matlab:

%Modified Secant Method

P = (0:0.1:20000);

l_k = 50;

Delta = 250;

eck = 0.4;

E = 200000;

lk = 50;

f_x = @(P) (Delta)./(1+((eck)*sec(0.5*sqrt(P/(E))*l_k)));

plot(P,f_x(P));

for iter = 1:3

    

    sol(iter) = P-(Delta*P*f_x(P))/(f_x(P+Delta*P)-f_x(P));

    Delta = P;

    P = sol(iter);

    

end

sol()

8.17

 8.30

8.46

Homework # 6

6.2

(a) Graphically I calculated the root to 3.56

Matlab: 

x = (3:0.001:4);

f_x = @(x) 2*x.^3 - 11.7*x.^2 + 17.7*x - 5;

plot(x,f_x(x));

(b) Using the Fixed-Point iteration method with up to three iterations I received this output:

Output:

3.5353    3.5514    3.5583

Matlab: 

%Fixed point iteration

x = 0;

f_x = @(x) (2*x.^3 - 11.7*x.^2 - 5)./-17.7;

sol = f_x(x)

for iter = 1:2

   sol = f_x(sol)  

end

(c) Using the Newton Raphson Method:

Output:

3.5672    3.5632    3.5632

Matlab:

%Newton Raphson Method
x = 3;

f_x = @(x) (2*x.^3 - 11.7*x.^2 + 17.7*x -5);

f_px = @(x) (6*x.^2 - 23.4*x +17.7);

sol(1) = x - (f_x(x))/(f_px(x));

for iter = 2:3
  x= sol(iter-1);
  sol(iter) = x - (f_x(x))/(f_px(x));

end

sol()

(d)Using the Secant method:

Output:

2.7463    2.9326    3.0570

Matlab:

  

%Secant Method

x = 3;

x_1 = 2;

f_x = @(x) (2*x.^3 - 11.7*x.^2 - 5)./-17.7;

for iter = 1:8

    sol(iter) = (f_x(x)*(x_1-x))/(f_x(x_1)-f_x(x));

    x_1 = x;

    x = sol(iter);

end

sol()

(e)Using the Modified Secant method:

Output:

3.1621    3.4121    3.6067

%Modified Secant Method

x = 3;

Delta = 4;

f_x = @(x) (2*x.^3 - 11.7*x.^2 - 5)./-17.7;

for iter = 1:3

    sol(iter) = x-(Delta*x*f_x(x))/(f_x(x+Delta*x)-f_x(x));

    Delta = x;

    x = sol(iter);

end

sol()

6.5

  

For this problem I just reused the matlab from problem 6.2, 4.2 Gave me a valid Result, however the input of 4.43 did not.

Output:

4.2 : 0.4746

4.43: -0.2280

Matlab:

%Newton Raphson Method

x = 4.2;

f_x = @(x) (-2 + 6*x -4*x.^2 +0.5*x.^3);

f_px = @(x) (1.5*x.^2 -8*x +6);

sol(1) = x - (f_x(x))/(f_px(x));

for iter = 2:8

  x= sol(iter-1);

  sol(iter) = x - (f_x(x))/(f_px(x));

end

sol()

 6.15

(b) Using this method(in my Matlab) I was able to execute a modified version of the Newton -Raphson method to calculate the roots of these equations:

Matlab:

%Newton Raphson Method for a non linear system

%u_0 =  x^2 +xy - 10= 0 

%v_0 =  y + 3xy^2 -57 = 0

x = 2;

y = 1.5;

for iter = 1 : 10

u_0 = (-x.^2 + x + 0.75)/y;

v_0 = (y+5*x*y)/x.^2;

du_0x = (1-2*x)/y;

du_0y = (x.^2 - x - 0.75)/y.^2;

dv_0x = -((5*x+2)*y)/x.^3;

dv_0y = (5*x +1)/x.^2;

deter = du_0x*dv_0y - du_0y*dv_0x;

xsol = x - ((u_0*dv_0y - v_0*du_0y)/(deter));

ysol = y - ((v_0*du_0x - u_0*dv_0x)/(deter));

sol1(iter) = xsol;

sol2(iter) = ysol;

x = xsol;

y = ysol;

end

sol1()

sol2()

6.18

For this problem I took the supplied equation and rewrote it as c = ((W-Qc)/kV)^2 and applied my modified secant method with c = 4 and delta = 0.5 and got an Output as follows:

Output:

7.6000    7.3750    7.2394

-0.0305   -0.0187

Matlab:

c = 4;

Delta = 0.5;

V = 1*10^6;

Q = 1*10^5;

W = 1*10^6;

k = 0.25;

f_c = @(c) ((W - Q*c)/k*V)^2;

for iter = 1:3

    sol(iter) = c-(Delta*c*f_c(c))/(f_c(c+Delta*c)-f_c(c));

    Delta = c;

    c = sol(iter);

end

sol()

err(1) = (sol(2)-sol(1))/sol(2)

err(2) = (sol(3)-sol(2))/sol(3)

6.28

Using the value of 16.15 will not give you the root between 15 and 20 because 16.15 is lower than the root, so it converges to the next root, 0.468. However if I change the initial guess to lets say 20, the root is found at : 18.8948

Matlab:

%Newton Raphson Method

x = 16.15;

f_x = @(x) (0.0074*x.^4 - 0.284*x.^3 + 3.355*x.^2 - 12.183*x + 5);

f_px = @(x) (0.0296*x.^3 - 0.852*x.^2 + 6.71*x - 12.183);

for iter = 1:9

  

  sol(iter) = x - (f_x(x))/(f_px(x));

  x = sol(iter);

end

sol()