Problem 4.25

Nasser Abbasi

 

Output

» nma_HW_4_25_newtn

 

  nma_HW_4_24_newtn - Program to solve a system of nonlinear equations

  using Newton's method.  Equations defined by function fnewt.

 

Enter the initial guess (row vector): [0 0.1981]

Enter number of step to iterate [N]:30

Enter name of your fnewt M file to be called:'nma_HW_4_24_fnewt'

k2 now =0.666667, step number=0 position now is x=0.0066564 m, y=0.191335 m

k2 now =1.33333, step number=0 position now is x=0.0126819 m, y=0.185281 m

k2 now =2, step number=0 position now is x=0.0179529 m, y=0.179878 m

k2 now =2.66667, step number=0 position now is x=0.0225787 m, y=0.175035 m

k2 now =3.33333, step number=0 position now is x=0.026657 m, y=0.170676 m

k2 now =4, step number=0 position now is x=0.0302691 m, y=0.166735 m

k2 now =4.66667, step number=0 position now is x=0.0334821 m, y=0.163161 m

k2 now =5.33333, step number=0 position now is x=0.0363522 m, y=0.159907 m

k2 now =6, step number=0 position now is x=0.0389261 m, y=0.156935 m

k2 now =6.66667, step number=0 position now is x=0.0412434 m, y=0.154214 m

k2 now =7.33333, step number=0 position now is x=0.0433372 m, y=0.151713 m

k2 now =8, step number=0 position now is x=0.0452357 m, y=0.14941 m

k2 now =8.66667, step number=0 position now is x=0.046963 m, y=0.147283 m

k2 now =9.33333, step number=0 position now is x=0.0485394 m, y=0.145313 m

k2 now =10, step number=0 position now is x=0.0499825 m, y=0.143486 m

k2 now =10.6667, step number=0 position now is x=0.0513074 m, y=0.141786 m

k2 now =11.3333, step number=0 position now is x=0.052527 m, y=0.140202 m

k2 now =12, step number=0 position now is x=0.0536528 m, y=0.138723 m

k2 now =12.6667, step number=0 position now is x=0.0546945 m, y=0.137338 m

k2 now =13.3333, step number=0 position now is x=0.0556607 m, y=0.136041 m

k2 now =14, step number=0 position now is x=0.0565589 m, y=0.134822 m

k2 now =14.6667, step number=0 position now is x=0.0573955 m, y=0.133676 m

k2 now =15.3333, step number=0 position now is x=0.0581766 m, y=0.132596 m

k2 now =16, step number=0 position now is x=0.058907 m, y=0.131576 m

k2 now =16.6667, step number=0 position now is x=0.0595915 m, y=0.130613 m

k2 now =17.3333, step number=0 position now is x=0.060234 m, y=0.129702 m

k2 now =18, step number=0 position now is x=0.0608381 m, y=0.128838 m

k2 now =18.6667, step number=0 position now is x=0.061407 m, y=0.128019 m

k2 now =19.3333, step number=0 position now is x=0.0619436 m, y=0.127241 m

k2 now =20, step number=0 position now is x=0.0624505 m, y=0.126501 m

After 30 iterations the final position is

    0.0625    0.1265

 

»

 


Source code

% nma_HW_4_25_newtn - Program to solve a system of nonlinear equations

% using Newton's method.  Equations defined by function fnewt.

clear all;  help nma_HW_4_25_newtn;  % Clear memory and print header

 

% Modified by Nasser Abbasi for problem 4.25

 

%* Set initial guess and parameters

x0      = input('Enter the initial guess (row vector): ');

 

x       = x0;  % Copy initial guess

xp(:,1) = x(:); % Record initial guess for plotting

 

nStep   = input('Enter number of step to iterate [N]:');

 

fnewt   = input('Enter name of your fnewt M file to be called:');

 

%* Loop over desired number of steps

 

%

% use continuation with k2 as the variable.

% we know that when k2=0, solution is x* = [0 0.1981]

% for each iteration i use

%

% k_i = k0 +  ( k_20 - k0 )  i/N

%

%

 

k0 = 0;

ka = 20;

 

   for iStep=1:nStep

 

      %

      % apply continuatin to find next value for k2

      %

      k = k0 +  (  ka - k0 ) * ( iStep / nStep );

   

      %* Evaluate function f and its Jacobian matrix D

      [f D] = feval(fnewt,x,k);      % fnewt returns value of f and D

      %* Find dx by Gaussian elimination

      dx = f/D;

      %* Update the estimate for the root 

      x = x - dx;              % Newton iteration for new x

      xp(:,iStep+1) = x(:); % Save current estimate for plotting 

 

      fprintf('k2 now =%g, step number=%d ',k,i);

      fprintf('position now is x=%g m, y=%g m\n',x(1),x(2));

   end

 

%* Print the final estimate for the root

fprintf('After %g iterations the final position is \n',nStep);

disp(x);

 

 

 


function [f,D] = nma_HW_4_24_fnewt(state,k2)

%  Function used by the N-variable Newton's method

%  Inputs

%    x         State vector [x y]

%    param     k2

%  Outputs

%    f     r.h.s. [dx/dt dy/dt]

%    D     Jacobian matrix, D(i,j) = df(j)/dx(i)

 

% Modified from original fnewtn to solve problem 4.24

% Nasser abbasi

 

% Evaluate f(i)

x= state(1);

y= state(2);

f(1) =  10*((x^2+y^2)^(1/2)-1/10)/(x^2+y^2)^(1/2)*x+...

        1/20*k2*(1/10*(100*x^2-20*x+1+100*y^2)^(1/2)-1/10)...

        /(100*x^2-20*x+1+100*y^2)^(1/2)*(200*x-20);

f(2) = 10*((x^2+y^2)^(1/2)-1/10)/(x^2+y^2)^(1/2)...

       *y+10*k2*(1/10*(100*x^2-20*x+1+100*y^2)^(1/2)-1/10)/...

       (100*x^2-20*x+1+100*y^2)^(1/2)*y-981/1000;

 

% Evaluate D(i,j)

D(1,1) = 10/(x^2+y^2)*x^2-10*((x^2+y^2)^(1/2)-...

         1/10)/(x^2+y^2)^(3/2)*x^2+10*((x^2+y^2)^(1/2)-...

         1/10)/(x^2+y^2)^(1/2)+1/400*k2/(100*x^2-20*x+1+100*y^2)*...

         (200*x-20)^2-1/40*k2*(1/10*(100*x^2-20*x+1+100*y^2)^(1/2)-...

         1/10)/(100*x^2-20*x+1+100*y^2)^(3/2)*(200*x-20)^2+...

         10*k2*(1/10*(100*x^2-20*x+1+100*y^2)^(1/2)-...

         1/10)/(100*x^2-20*x+1+100*y^2)^(1/2);

 

D(1,2) = 10/(x^2+y^2)*y*x-10*((x^2+y^2)^(1/2)-...

         1/10)/(x^2+y^2)^(3/2)*x*y+1/2*k2/(100*x^2-20*x+1+100*y^2)...

         *y*(200*x-20)-5*k2*(1/10*(100*x^2-20*x+1+100*y^2)^(1/2)...

         -1/10)/(100*x^2-20*x+1+100*y^2)^(3/2)*(200*x-20)*y;

 

D(2,1) = 10/(x^2+y^2)*y*x-10*((x^2+y^2)^(1/2)-1/10)/(x^2+y^2)^(3/2)*x*y+...

         1/2*k2/(100*x^2-20*x+1+100*y^2)*y*(200*x-20)-...

         5*k2*(1/10*(100*x^2-20*x+1+100*y^2)^(1/2)-1/10)/(100*x^2-20*x+...

         1+100*y^2)^(3/2)*(200*x-20)*y;

 

D(2,2) = 10/(x^2+y^2)*y^2-10*((x^2+y^2)^(1/2)-1/10)/(x^2+y^2)^(3/2)*y^2+...

         10*((x^2+y^2)^(1/2)-1/10)/(x^2+y^2)^(1/2)+100*k2/(100*x^2-20*x+1+...

         100*y^2)*y^2-1000*k2*(1/10*(100*x^2-20*x+1+100*y^2)^(1/2)-...

         1/10)/(100*x^2-20*x+1+100*y^2)^(3/2)*y^2+10*k2*(1/10*(100*x^2-20*x...

         +1+100*y^2)^(1/2)-1/10)/(100*x^2-20*x+1+100*y^2)^(1/2);

 

 

return;