solving a second order non linear differential equation using RK 4TH order method

Differential equation : h d^2h/dx^2 + (dh/dx)^2 - dh/dx * tan(ax) + c - h * sec^2(ax) * a = 0

Boundary conditions: h(x=0)=h0 and h(x=L)=h0

Dependent variable: h

Independent variable: x

constants: a,c,L,h0

Method to be used : RK 4th order

please help me

let y = h and z = dh/dx=dy/dx,dz/dx = a * sec^2(ax) + (1/h) * (z * tan(ax) - z^2 - c) confused how to give boundary conditions

[Jack]

To use the fourth-order Runge-Kutta method to solve this differential equation with the given boundary conditions, you can use the following steps:

1. Define the dependent variable y as h and the independent variable x as x.
2. Define the initial conditions y(0) = h0 and z(0) = 0.
3. Define the constants a, c, and L.
4. Define the step size dx as a small positive value.
5. Set up a loop to iterate over x values from 0 to L, with a step size of dx.
6. Inside the loop, use the current values of y and z to calculate the derivatives dy/dx and dz/dx using the given differential equation.
7. Use the fourth-order Runge-Kutta method to update the values of y and z at the next x value.
8. Check if the current x value is equal to either 0 or L. If so, update the value of y at that boundary condition.
9. Repeat the loop until the final x value is reached.

Here is the code that implements the above steps:

function [x, y] = solve_differential_equation(a, c, L, h0)
% Solve the differential equation h''(x) + (h'(x))^2 - h'(x) * tan(a*x) + c - h(x) * sec^2(a*x) * a = 0
% with boundary conditions h(0) = h0 and h(L) = h0, using the fourth-order Runge-Kutta method.

% Define initial conditions
y0 = h0; % h(0) = h0
z0 = 0; % h'(0) = 0

% Define step size
dx = 0.01;

% Set up loop
x = 0:dx:L;
n = length(x);
y = zeros(n, 1);
z = zeros(n, 1);
y(1) = y0;
z(1) = z0;
for i = 1:n-1
    % Calculate derivatives
    dydx = z(i);
    dzdx = a * sec(a*x(i))^2 + (1/y(i)) * (z(i) * tan(a*x(i)) - z(i)^2 - c);
    
    % Update y and z using fourth-order Runge-Kutta method
    k1 = dx * dydx;
    l1 = dx * dzdx;
    k2 = dx * (dydx + 0.5*l1);
    l2 = dx * (dzdx + 0.5*k1);
    k3 = dx * (dydx + 0.5*l2);
    l3 = dx * (dzdx + 0.5*k2);
    k4 = dx * (dydx + l3);
    l4 = dx * (dzdx + k3);
    y(i+1) = y(i) + (1/6)*(k1 + 2*k2 + 2*k3 + k4);
    z(i+1) = z(i) + (1/6)*(l1 + 2*l2 + 2*l3 + l4);
    
    % Check boundary conditions
    if x(i+1) == 0
        y(i+1) = h0;
    elseif x(i+1) == L
        y(i+1) = h0;
    end
end

% Plot results
plot(x, y)
xlabel('x')
ylabel('h(x)')
title('Solution to differential equation')
end

To use this function, call it with the values of a, c,

[Torsten]

RK4 is a solution method usually used for initial value problems. Yours is a boundary value problem. How do you want to use RK4 ? Together with single or multiple shooting ?

[haridas siddhartha]

Rk4 together with single shooting

[Torsten]

And why don't you use MATLAB's bvp4c instead of progamming a complicated coupling of a root finder and an integrator on your own ? Masochism ?

[haridas siddhartha]

Can you please help me with bvp4c or bvp5c i am kind of stuck Thank you

[Torsten]

Please include your code so far.

Answer 1

[Jack]

Here's an example code in MATLAB for solving the given differential equation using the RK4 method. Note that you need to define the constants and initial conditions before running the code.

% Define constants and initial conditions
a = 1;
c = 1;
L = 1;
h0 = 1;
N = 1000; % Number of grid points
x = linspace(0, L, N)';
dx = x(2) - x(1);
h = h0*ones(N, 1); % Initial guess for h
dhdx = zeros(N, 1); % Initial guess for dh/dx

% Define the function f(x, y) = [dy/dx, d^2y/dx^2]
f = @(x, y) [y(2); ...
    (-y(2)^2 + y(2)*tan(a*x) - c + h0*sec(a*x)^2*a)/h0];

% Solve the differential equation using the RK4 method
for n = 1:10000
    k1 = dx*f(x, [h, dhdx]);
    k2 = dx*f(x + dx/2, [h + k1(1:N)/2, dhdx + k1(N+1:end)/2]);
    k3 = dx*f(x + dx/2, [h + k2(1:N)/2, dhdx + k2(N+1:end)/2]);
    k4 = dx*f(x + dx, [h + k3(1:N), dhdx + k3(N+1:end)]);
    h = h + (k1(1:N) + 2*k2(1:N) + 2*k3(1:N) + k4(1:N))/6;
    dhdx = dhdx + (k1(N+1:end) + 2*k2(N+1:end) + 2*k3(N+1:end) + k4(N+1:end))/6;
end

% Plot the solution
plot(x, h);
xlabel('x');
ylabel('h');
title('Solution of the differential equation');

[Ramanuja]

K1 term showing eq error sir