Solve 3в order ODE

Hello everyone!

Im trying to solve following equation f(x):

(1e-7+5.1471e-04+1.164e2*f')*f''' +1.164e2*f'''^2-1.2e-06*f'-12/5*(f''*(1+f)+f'^2)=0

with

f(0)= -1

f(xmax)=0

f'(0)=0

I have no clue how I can rearange it since there are f'^2 and f''*f terms.

As far as I know the ODE system should look like:

df(1)=f'...

df(2)=f''...

df(3)=f'''...

df(4)=f''''...

but I dont know how I can make it work.

Any help would be really helpfull...

[John D'Errico]

I don't see any 4th order term in there. Unless, perhaps that is what you intend by f''*.

So assuming that to be true, a 4th order ODE must have exactly 4 initial conditions. You only seem to be providing 3 of them. And one of them is a condition at the far end, so a boundary condition. It also means you will need to use some scheme other then the initial value solvers normally used. A shooting method is one common approach, as is collocation, but even those must fail unless you provide a 4th condition. (Assuming this really is a 4th order ODE.)

I'm not even sure that what we see is the real problem, since I see things like:

1e-7+5.1471e-04

Is there a reason why you would not add those two constants? My question is if there are some other terms attached to them.

Anyway, since this is a fairly nonlinear ODE, you will not easily arrange it into a simple form as you think you want to do. I would guess a collocation scheme might be most appropriate.

[Matdoomer]

sorry, its a 3d oder problem, not 4th!

Answer 1

[Ameer Hamza]

To solve an n-th order ODE functions with MATLAB, you need to express your equation in the form

Since your equations are quadratic in term of so you will have two different solutions for your equation. See the following code. It uses quadratic formula to write in term of and

% (1e-7+5.1471e-04+1.164e2*f')*f''' +1.164e2*f'''^2-1.2e-06*f'-12/5*(f''*(1+f)+f'^2)=0

xmin = 0;
xmax = 1;
x = linspace(xmin, xmax, 20);
init = bvpinit(x, [0;0;0]);
sol1 = bvp4c(@odeFun1, @bcFun, init);
sol2 = bvp4c(@odeFun2, @bcFun, init);

figure;
plot(sol1.x, real(sol1.y))
title('Solution 1')
legend({'f', 'f''', 'f'''''});

figure;
plot(sol2.x, real(sol2.y))
title('Solution 2')
legend({'f', 'f''', 'f'''''});

function dfdx = odeFun1(x, f)
    dfdx(1) = f(2);
    dfdx(2) = f(3);
    
    % write as quadratic equation of form a*x^2+b*x+c=0
    a = 1.164e2;
    b = (1e-7+5.1471e-04+1.164e2*f(2));
    
    c = -1.2e-06*f(2)-12/5*(f(3)*(1+f(1))+f(2)^2);
    dfdx(3) = (-b + sqrt(b^2-4*a*c))/(2*a);
end

function dfdx = odeFun2(x, f)
    dfdx(1) = f(2);
    dfdx(2) = f(3);
    
    % write as quadratic equation of form a*x^2+b*x+c=0
    a = 1.164e2;
    b = (1e-7+5.1471e-04+1.164e2*f(2));
    
    c = -1.2e-06*f(2)-12/5*(f(3)*(1+f(1))+f(2)^2);
    dfdx(3) = (-b - sqrt(b^2-4*a*c))/(2*a);
end

function res = bcFun(xa, xb)
    res = [xa(1)+1; % f(0)= -1
           xb(1);   % f(xmax)=0
           xa(2)];  % f'(0)=0        
end