RK$ with shooting algorithm for two boundary value problem

1 次查看（过去 30 天）
Debayan Dasgupta
Debayan Dasgupta 2022-5-11
回答： Nipun 2024-6-3
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.
  4 个评论
显示 2更早的评论
Debayan Dasgupta
Debayan Dasgupta 2022-5-12
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
Torsten 2022-5-12
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Nipun
Nipun 2024-6-3
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

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by