Numerical Analysis. Please, help with MATLAB.

Question

Use the initial approximation $(p_0,q_0)=(-0.3,-1.3)$, and compute the next three approximations to the fixed point using

a) Fixed-point iteration and equations : $p_{k+1} = g_{1}(p_k,q_k)$ and $q_{k+1} = g_2(p_k,q_k)$

b) Seidel iteration using equations $p_{k+1}=g_1(p_k,q_k)$ and $q_{k+1} = g_2(p_{k+1},q_k)$

$$x=g_1(x,y) = \dfrac{y-x^3+3x^2+3x}{7} \text{ (cubic)}$$

$$y = g_2(x,y)=\dfrac{y^2+2y-x-2}{2} \text{ (parabola).}$$

So, the direct computations gave me: a)Fixed point iteration:

$ (p_0,q_0) = (-0.3,-1.3)$

$(p_1,q_1) = (-0.2684,-1.3175)$

$(p_2,q_2) = (-0.2694,-1.3161)$

$(p_3,q_3) = (-0.2696, -1.3153)$

b)Seidel iteration:

$ (p_0,q_0) = (-0.3,-1.3)$

$(p_1,q_1) = (-0.2719,-1.3191)$

$(p_2,q_2) = (-0.2704,-1.3139)$

$(p_3,q_3) = (-0.2694, -1.3160)$

I am wondering if there is a way to program it in MATLAB, so I could get more elegant solution.

Answer 1

I'm no pro at MATLAB, but I know some basics. You should be able to do it like this. First define a function, here I did it for you seidels iteration, because you didnt provide equations for the other one.

function [ pk qk ]= seidel_iteration(x,y)

pk=(y-x^3+3*x^2+3*x)/7;

qk=(y^2+2*y-x-2)/2;

endfunction

Then you can call the function to print out the values you want

p1=-0.3

q1=-1.3

numberIterations= 10 % or whatever number of iterations you want

for i=1:1:numberIterations

disp([ p1 q1])

[ p1 q1 ]= seidel_iteration(p1,q1);

endfor

I am truly sorry about the formatting, it is horrible. I dont know how to write code nicely into Math Stack

Answer 2

This is a slightly modified version, working with Matlab 2010. It is based on the code of @Sebastian. All credit goes to him.

function fp_iter
clc        % I like a clean screen
p1=-0.3;
q1=-1.3;
tol=1e-10; % tolerance for convergence
i=0;


while true
    i=i+1; % count number of iteration

    str=sprintf('At iteration %d we have:   \t q1=%0.10f and p1=%0.10f',i,p1,q1);
    disp(str);

    [ newp1 newq1 ]= seidel_iteration(p1,q1); % perform iteration
    if abs(newp1-p1)+abs(newq1-q1)<tol        % check if converged
        disp('We reached the required tolerance');
        break
    end

    res=sprintf('Difference compared to previous iteration is %0.10f\n',...
        abs(newp1-p1)+abs(newq1-q1));
    disp(res);

    p1 = newp1; % update parameter p1
    q1 = newq1; % update parameter q1
end

end
% if you put an "end" after each function, multiple function can be in one file.

function [ pk qk ]= seidel_iteration(x,y)

pk=(y-x^3+3*x^2+3*x)/7;
qk=(y^2+2*y-x-2)/2;

end