Shooting Method with Secant Method

Hi all,

I am working on a problem to solve a second order differential equation. I am using the numerical method 'the shooting method' and I need to perform iterations to find the initial value of the slope. To do this I am using the secant method. My intial values to start the shooting method for z are -4 and 4. From this is get the corresponding y values of -65.025 (3dp) and 106.062 (3dp) respectively. I am aiming to find the value of z when y = 0. Thus, I proceed to the secant method to find the value of 'z' when y = 0;

The formula for secant method is:

z2 = z1 - (y(z1) - 0)*(z1 - z2)/(y(z1 - y(z0))

I want to make a loop where this updates: Hence:

z0 = z1

z1 = z2

y(z0) = y(z1)

y(z1) = y(z2)

So, the next iteration is:

z3 = z2 - (y(z2) - 0)*(z2 - z21/(y(z2 - y(z2))

The thing is I now need to compute y(z2). I can do this, but it would meaning copying and pasting another block of code. Thus, my question is how to make the algorithm in MATLAB so that I can perform the desired number of iterations I want within the for loop? Rather than copying and pasting the code below over and over again until I reach a satisfactory point of convergence?

Here is my code:

clear, clc, close all

format longG

% Structural Properties
E = 2.1E+08;
Ic = 22530/100^4;
Ib = 19460/100^4;
Ac = 168/100^2;
h = 0.5;
r = 1;
Gi = (1*E*Ib/10);
Ci = 1/2*(E*Ic/h);
g = Gi/Ci;
K = (6*g + 1 + r);
V = [0 0 0 0 0 35 0 0 0 0 0 25 0 0 0 0 0 15 0 0 0 0 0 5 0];

%Initial Conditions 1
y = 0
z = -4
             for i = 1:24
                    if i == 6 || i == 12 || i == 18
                    F(i) = h*(V(i) + V(i+6))/(4*Ci);
                    z(i+1) = z(i) + ((F(i) - K*z(i) + y(i))/(-z(i)))*h;
                    y(i+1) = y(i) + z(i)*h;
             else F(i) = 0;
                    z(i+1) = z(i) + ((F(i) - K*z(i) + y(i))/(-z(i)))*h;
                    y(i+1) = y(i) + z(i)*h;
             end
             if i == 24
                   F(i) = h*(V(i))/(4*Ci);
                   z(i+1) = z(i) + ((F(i) - K*z(i) + y(i))/(-z(i)))*h;
                   y(i+1) = y(i) + z(i)*h
             end
             q0 = y(end)
             p0 = z(1)
             end
q0 = y(end)
p0 = z(1)

%Initial Conditions 2
y = 0
z = 4
             for i = 1:24
                    if i == 6 || i == 12 || i == 18
                    F(i) = h*(V(i) + V(i+6))/(4*Ci);
                    z(i+1) = z(i) + ((F(i) - K*z(i) + y(i))/(-z(i)))*h;
                    y(i+1) = y(i) + z(i)*h;
             else F(i) = 0;
                    z(i+1) = z(i) + ((F(i) - K*z(i) + y(i))/(-z(i)))*h;
                    y(i+1) = y(i) + z(i)*h;
             end
             if i == 24
                   F(i) = h*(V(i))/(4*Ci);
                   z(i+1) = z(i) + ((F(i) - K*z(i) + y(i))/(-z(i)))*h;
                   y(i+1) = y(i) + z(i)*h
             end
             end
q1 = y(end)
p1 = z(1)

% Secant Method

p2 = p1 - (q1 - 0)*(p1 - p0)/(q1 -q0)

I have left the secant method formula outside the loop, because I am unsure how to include it in the loop to perform the iterations I previously explained.

I understand that this is quite a long question. So thank you in advance for any help. However, you guys on here are great, so I hoping its not too tedious for you.

Many thanks,

Scott

Answer 1

[surya venu]

Hi,

To integrate the secant method into your MATLAB code and perform the iterations automatically, you can wrap the entire process in a loop. This way, you don't have to manually repeat the shooting method calculations for new guesses of "z". Here's how you can structure your code to do this:

1. Define a maximum number of iterations and a tolerance for convergence. The tolerance will help you determine when the value of "y" is close enough to 0 to stop the iterations.
2. Use a loop to repeatedly apply the shooting method with new guesses for "z" based on the secant method until either the maximum number of iterations is reached or the convergence criterion is met.
3. Update the values of z0, z1, y(z0), and y(z1) at the end of each iteration as per the secant method.

Here's a modified version of your code incorporating these suggestions:

clear, clc, close all
format longG

% Structural Properties
E = 2.1E+08;
Ic = 22530/100^4;
Ib = 19460/100^4;
Ac = 168/100^2;
h = 0.5;
r = 1;
Gi = (1*E*Ib/10);
Ci = 1/2*(E*Ic/h);
g = Gi/Ci;
K = (6*g + 1 + r);
V = [0 0 0 0 0 35 0 0 0 0 0 25 0 0 0 0 0 15 0 0 0 0 0 5 0];

% Initial guesses and values
z0 = -4;
z1 = 4;
maxIters = 100; % Maximum number of iterations
tol = 1e-5; % Convergence tolerance
iter = 0;

while iter < maxIters
    y = zeros(1,25); % Reset y for each shooting method iteration
    z = zeros(1,25); % Reset z for each shooting method iteration
    z(1) = z0; % Set initial condition for z
    
    % Shooting method for z0
    for i = 1:24
        if i == 6 || i == 12 || i == 18
            F(i) = h*(V(i) + V(i+6))/(4*Ci);
        else
            F(i) = 0;
        end
        z(i+1) = z(i) + ((F(i) - K*z(i) + y(i))/(-z(i)))*h;
        y(i+1) = y(i) + z(i)*h;
    end
    q0 = y(end);

    % Reset y and z for the next initial condition
    y = zeros(1,25);
    z = zeros(1,25);
    z(1) = z1; % Set initial condition for z
    
    % Shooting method for z1
    for i = 1:24
        if i == 6 || i == 12 || i == 18
            F(i) = h*(V(i) + V(i+6))/(4*Ci);
        else
            F(i) = 0;
        end
        z(i+1) = z(i) + ((F(i) - K*z(i) + y(i))/(-z(i)))*h;
        y(i+1) = y(i) + z(i)*h;
    end
    q1 = y(end);
    
    % Check for convergence
    if abs(q1) < tol
        fprintf('Converged to y = 0 with z = %f after %d iterations\n', z1, iter);
        break;
    end
    
    % Update for next iteration using secant method
    p2 = z1 - (q1 - 0)*(z1 - z0)/(q1 - q0);
    
    % Prepare for next iteration
    z0 = z1;
    z1 = p2;
    q0 = q1;
    
    iter = iter + 1;
end

Converged to y = 0 with z = 5.407086 after 18 iterations

if iter == maxIters
    fprintf('Max iterations reached without convergence.\n');
end

Hope it helps.

[Scott Banks]

Thankyou, that is fantastic!

On a further note, however, if I change the Ci variable that also has 'h' in it - if I change this to 3 (which is what it should be), the solution does not converge, even after thousands of iterations.

Do you know why this may be?