Solving 4th order ode using ode45

How to solve the following 4th order ode using ode45 solver

[madhan ravi]

Are you sure ode45() is the one to use? Seems like a BVP problem.

[Bjorn Gustavsson]

...and aren't you missing one boundary-value?

@madhan - perhaps just as "easy" to use shooting method to find a solution that trails of towads flat at infinity?

[Elayarani M]

Thank you. I will try with that method.

Answer 1

[Alan Stevens]

Here's some coding that basically solves the equation. I've no idea what the value of k should really be, but the constants chosen give a consistent result. The choice of f'''(0) is based on the original equation with the other x=0 values plugged in; where f''(0) is a chosen to give a seemingly reasonable result!

k = -0.002;
xspan = [0 100];
d2fdx20 = -1;
F0 = [0 1 d2fdx20 (1-k*(d2fdx20^2))/(1-2*k)];  

[x, F] = ode45(@rates, xspan, F0, [], k);

f = F(:,1);
dfdx = F(:,2);
plot(x, f, x, dfdx),grid
xlabel('x'), ylabel('f and dfdx')
legend('f','dfdx')

function dFdx = rates(x,F,k)
        
         f = F(1);
         dfdx = F(2);
         d2fdx2 = F(3);
         d3fdx3 = F(4);
         if x==0
            d4fdx4 = 0;
         else
            d4fdx4 = (d3fdx3 +f.*d2fdx2 - dfdx.^2 - 2*k*dfdx.*d3fdx3 + k*d2fdx2.^2)./(k*f);
         end
         
         dFdx = [dfdx; d2fdx2; d3fdx3; d4fdx4];
         
end

[Alan Stevens]

In the above f'(x) actually goes to a small negative value, which gradually drags f down. If we assume the small negative value is a result of numerical inaccuracies, then by replacing

dfdx = F(2);

by

dfdx = max(F(2),0);

in function "rates", we get f tending to a constant.

[Elayarani M]

Thank you. I have added a negative sign in d4fdx4 expression and taken positive value for k. It works well.

[Elayarani M]

Hi Alan. when f(0)=1, how to take F0?

[Bjorn Gustavsson]

If f is a function, then just like that, if f is the output from your ODE, then it is an array where each element matches the function-values (f and its derivatives) at the values of x. Then you get "f(x=0)" as: F(1,1) if x(1) == 0.

[Alan Stevens]

Set the first element of F0 to 1.

[Elayarani M]

Yes, first element of F0 is 1. What will be the fourth element?

[Alan Stevens]

Assuming that f''''(0) = 0, then use your ODE with f(0) = 1 and rearrange it to get f'''(0) in terms of f(0), f'(0), f''(0) and k. You will still have to make an assumption about f''(0) though (I'd be inclined to keep it at -1 initially).

[Elayarani M]

Is it correct to assume fourth derivative of f to zero at x=0. Will it not affect the result?

[Alan Stevens]

The expression for the fourth derivative is

d4fdx4 = (d3fdx3 +f.*d2fdx2 - dfdx.^2 - 2*k*dfdx.*d3fdx3 + k*d2fdx2.^2)./(k*f);

This has f in the denominator. When f(0) is zero, the only way of avoiding a singularity at x = 0 is to assume the derivative itself is zero. This isn't necessarily the case when f(0) = 1, so you could remove the part of the " if" statement that sets it to zero in this case.

[Elayarani M]

Yes, I have removed that part for f(0)=1. Should i put f""(0)=0 for finding f'''(0) expression in F0 array? Is it correct to put that?

[Alan Stevens]

Yes if it leads to the numerator of d4fdx4 being zero; then it would be consistent at least.

[Elayarani M]

Thank you Alan Stevens. I will try with that.

[Elayarani M]

For the case f(0)=1

If i put f''''(0)=0 for finding f'''(0) expression, am not getting the correct value. I tried with two assumptions but its not working correctly. I dont know how to take f'''(0) expression in F0 array. Can anyone help?

[Alan Stevens]

What's the value of k you are using?

[Bjorn Gustavsson]

For a first-order ODE you need to specify one boundary-condition, for a second-order ODE you need to specify 2 boundary conditions (e.g. with an equation of motion you need to specify the initial position and velocity). For an n-th order ODE you have to specify n boundary conditions. You have a 4th-order ODE with 3 specified boundary-conditions (two initial, for f(0) and f'(0) and one end-condition for f'(inf)). One way to handle this is to view this problem as a parameter-estimation-problem where you have to estimate the initial conditions for f'' and f''' that satisfies your end-condition, this method is known as the shooting method. If you use this you have to make a search for the values of f''(0) and f'''(0) that gives you a value of f'(inf) that fits.

[Elayarani M]

k=-1; d4fdx4 = -(d3fdx3 +f.*d2fdx2 - dfdx.^2 - 2*k*dfdx.*d3fdx3 + k*d2fdx2.^2)./(k*f); When f(0)=0, f''''(0) term vanishes in f'''(0) expression but when f(0)=1,f''''(0) term exists. I dont know how to handle that.

[Elayarani M]

Hi Bjorn Gustavsson. Using shooting method, when i make two initial guesses for the unknowns f''(0) and f'''(0) and make search for the values of f''(0) and f'''(0) that satisfies f'(inf)=0, it gives some value but not the correct value.

[Bjorn Gustavsson]

Do you get a solution that satisfies the ODE and your given initial conditions (f(0)=0, and f'(0)=1), and the end-condition?

If you need to get to infinity in a more proper sense you might have to make a variable-substitution, something like

t = tan(theta)

and then transform your ODE too.

[Elayarani M]

Yes, I get some solution that satisfies the given conditions. I will try with variable-substitution.

[Bjorn Gustavsson]

In the case that you get a solution that matches your boundary-conditions and satisfies your ODE why is that not good enough?

[Elayarani M]

The solution was not correct when compared with analytic solution.

[Bjorn Gustavsson]

Did the numerical solution differ by much? If not then perhaps only numerical deviations? Since you have a non-linear ODE there might be many solutions (right?), have you gotten all analytically? If not then the numerical solution might be one of the other valid solutions.

[Elayarani M]

Thank you Bjorn Gustavsson. First I will try to get all analytic solution and then cross check with the numerical solution.