Numerical Computation Matlab

Question

I am trying to solve the fifty problems $f_i(\theta) = 0, \ i = 1, ..., 50,$ where

$$f_i (\theta) = \theta - \sin(\theta) - \frac{i}{100}$$

I was told I must solve these in ascending order of $i$ and with a maximum of 200 evaluations of $f$ or its derivatives.

I tried to solve this in Matlab using Newton's method, with the following iteration formula:

$$\theta_{n+1} = \theta_n - \frac{f(\theta_n)}{f'(\theta_n)} = \theta_n - \frac{\theta_n- \sin(\theta_n) - \frac{i}{100}}{1- cos(\theta_n)}$$

Here's my basic Matlab code for one iteration applied to each problem:

x=0.01;
for i=1:1:50
    x=x - ((x-sin(x)-(i/100))/(1-(cos(x))))
end

However I think the question asks for 200 iterations per problem, so I modified the code to:

for i=1:1:50
    n_itn=200;
    x(1)=0.01;
    for k = 1:n_itn
    x(k+1)=x(k) - ((x(k)-sin(x(k))-(i/100))/(1-(cos(x(k)))))
    end
end

But I get an extremely long answer, here's a part of it:

  Columns 196 through 198

  -0.000004145710670  -0.000002034384641  -0.000000336803024

  Columns 199 through 201

  -0.000000117100873   0.000000102162712  -0.000150151084928

So, I am not sure which part of this is my solution. So, what is wrong with my code? Am I even on the right track?

Any help would be greatly appreciated.

Answer 1

I think the aim of this task is to teach you the wonders of the homotopy continuation method. Thus the goal is to solve $$ x-\sin x =\frac12 $$ which then gets subdivided into 50 steps raising the "watermark" of the right side from $0$ to $\frac12$, where the initial problem $$ x-\sin x=0 $$ has the easy (and unique) solution $x_0=0$.

The idea is that if $x_k$ is a good solution to $$ x-\sin x=\frac{k}{100} $$ then the error $f_{k+1}(x_k)=-\frac1{100}$ leads after two Newton iterations to an error of magnitude $(0.01)^{2^2}=10^{-8}$.

Two Newton iterations because 200 evaluations of $f$ and $f'$ for 50 problems leads to an average of 4 evaluations equal to 2 Newton steps.

Since the solution for $k=0$ is singular, one can not start the iteration for $x_1$ with $x=x_0=0$. Better approximate $x-\sin x\simeq\frac16x^3$ and thus start the iteration with $x=\sqrt[3]{\frac3{50}}=0.3914867641169…$ for $k=1$ and similar for further small $k$. Then use linear extrapolation to get good initial values. (Basic homotopy uses constant continuation.)

---

x = zeros(1,50);
for k=1:1:50
    if k<10 then
        x(k) = (3*k/50)^(1/3);
    else
        x(k) = 2*x(k-1)-x(k-2);
    end  
    for m = 1:2
        x(k) = x(k) - ( x(k)-sin(x(k))-k/100 ) / ( 1-cos(x(k))  );
    end
end

leading to results (partial list

k           x(k)          x(k) - sin(x(k))
-----------------------------------------
 1.    0.3924933889707    0.0100000000013
 2.    0.4952635544934    0.0200000000157
 3.    0.5676633468267    0.0300000000686
 4.    0.6255141630440    0.0400000001950
 5.    0.6745314424308    0.0500000004379
 6.    0.7175125729004    0.0600000008470
 7.    0.7560609279733    0.0700000014780
 8.    0.7911911702720    0.0800000023919
 9.    0.8235909073405    0.0900000036545
10.    0.8537501566741    0.1000000000114
15.    0.9811216100757    0.1500000000005
20.    1.0836918803146    0.2000000000001
25.    1.1712296525017    0.2500000000000
30.    1.2485154675427    0.3
35.    1.3182890721987    0.35
40.    1.382284133718     0.4
45.    1.4416754785973    0.45
50.    1.4973003890959    0.5

Answer 2

You should be iterating over the same x(i) therefore the index of x should be i not k+1. Hence:

x(i)=0.01;

and

x(i)=x(i) - ((x(i)-sin(x(i))-(i/100))/(1-(cos(x(i)))));

n_itn=200;
x = 0.01*ones(1,50);
for i=1:1:50    
    for k = 1:n_itn
        x(i)=x(i) - ((x(i)-sin(x(i))-(i/100))/(1-(cos(x(i)))));
    end
end