RK$ with shooting algorithm for two boundary value problem

I have a fourth order ODE of the form y''''=V*((x+U/V)*y'-y)/y^3-(3*y^2*y'''*y')/y^3; HEre U and V and theta are constants with known values, U=U=0.221951; V=0.04263; theta=0.2826;

For solving this, I need four initial conditions. But I have only two initial conditions and rest two are boundary condiitions. These are y(0)=1, y'(0)=0; y''(8)=0 and y'(8)=theta. Hece, I have proceeded with shooting algorithm along with RK4 to guess the intial conditions for y''(0) and y'''(0). But I am not able to satisfy both the conditions y''(8)=0 and y'(8)=theta at the same time. Any help will be highly appreciated.

[Jan]

How did you try it? What abpout BVP4C? Which code did you use? What does "I'm not able" mean?

[Torsten]

If you are not able to implement the shooting method for two unknown boundary values, you should use BVP4C, as Jan suggested.

But it's a simple Newton-update:

Estimate values for a = y''(0) and b = y'''(0).

Solve your ODE and get values fa = y'(8) and fb = y''(8).

Now solve your ODE for a+h = y''(0) and b = y'''(0) to get values faa = y'(8) and fba = y''(8).

Now solve your ODE for a = y''(0) and b+h = y'''(0) to get values fab = y'(8) and fbb = y''(8).

Now update a and b according to the usual Newton step

[a';b'] = [a;b] - [(faa-fa)/h, (fab-fa)/h; (fba-fb)/h,(fbb-fb)/h]^(-1) * [fa-theta;fb];

set a = a', b= b' and continue to do the same until

norm ([a-a';b-b']) < eps

and

norm([fa-theta;fb]) < eps

[Debayan Dasgupta]

Thank you so much for the suggestion. However, in Newton step, in the LHS a' or a has only one element i.e one value for each iteration, whereas the RHS has two elements in the form of (faa-fa)/h, (fab-fa)/h. This will create an error. I am missing something here?

[Torsten]

You mean the update

[a';b'] = [a;b] - [(faa-fa)/h, (fab-fa)/h; (fba-fb)/h,(fbb-fb)/h]^(-1) * [fa-theta;fb];

?

[a';b'] : 
2x1 vector

[a;b] : 
2x1 vector

[(faa-fa)/h, (fab-fa)/h; (fba-fb)/h,(fbb-fb)/h]^(-1) * [fa-theta;fb] : 
2x2 -matrix, multiplied by 2x1-vector = 2x1 vector.

All dimensions are compatible.

Answer 1

[Nipun]

Hi Debayan,

I understand that you are trying to solve the the fourth-order ODE using the shooting method and fourth order Range-Kutta method in MATLAB.

I assume that the initial conditions mentioned are correct and form the basis of the implementation. I recommend the following implementation to get the desired results:

function main
    % Constants
    U = 0.221951;
    V = 0.04263;
    theta = 0.2826;

    % Initial conditions
    y0 = [1; 0];  % y(0) = 1, y'(0) = 0
    x0 = 0;
    x_end = 8;

    % Initial guesses for y''(0) and y'''(0)
    guess1 = 0;
    guess2 = 0;
    tol = 1e-6;
    max_iter = 100;
    
    % Shooting method
    for iter = 1:max_iter
        % Solve ODE with current guesses
        [x, Y] = solve_ode(x0, x_end, [y0; guess1; guess2], U, V);

        % Calculate the boundary condition residuals
        res1 = Y(end, 3); % y''(8) should be 0
        res2 = Y(end, 2) - theta; % y'(8) should be theta

        % Check convergence
        if abs(res1) < tol && abs(res2) < tol
            break;
        end

        % Adjust guesses
        guess1 = guess1 - 0.1 * res1;
        guess2 = guess2 - 0.1 * res2;
    end

    if iter == max_iter
        warning('Max iterations reached without convergence');
    end

    % Plot the solution
    plot(x, Y(:, 1));
    xlabel('x');
    ylabel('y');
    title('Solution of the ODE');
end

function [x, Y] = solve_ode(x0, x_end, init_conditions, U, V)
    h = 0.01;
    x = x0:h:x_end;
    Y = zeros(length(x), 4);
    Y(1, :) = init_conditions';

    for i = 1:length(x)-1
        k1 = h * odefunc(x(i), Y(i, :), U, V);
        k2 = h * odefunc(x(i) + h/2, Y(i, :) + k1/2, U, V);
        k3 = h * odefunc(x(i) + h/2, Y(i, :) + k2/2, U, V);
        k4 = h * odefunc(x(i) + h, Y(i, :) + k3, U, V);
        Y(i+1, :) = Y(i, :) + (k1 + 2*k2 + 2*k3 + k4) / 6;
    end
end

function dYdx = odefunc(x, Y, U, V)
    y = Y(1);
    y1 = Y(2);
    y2 = Y(3);
    y3 = Y(4);
    y4 = V * ((x + U/V) * y1 - y) / y^3 - (3 * y^2 * y3 * y1) / y^3;
    dYdx = [y1; y2; y3; y4];
end

Explanation:

1. Main function:

• Defines constants, initial conditions, and initial guesses.
• Uses the shooting method to adjust y''(0) and y'''(0) until the boundary conditions are met.

2. solve_ode function:

• Uses RK4 to solve the ODE system from x0 to x_end.

3. odefunc function:

• Defines the system of first-order ODEs.

Hope this helps.

Regards,

Nipun