help solving differential equations

hello, i tring to solve this equation

m*f'^2 -0.5(m+1)*f*f'' = m+f'''

m is a variable i need to set

my initial values are:

f(0) = f'(0) = 0
f ' (inf) =1

this is my code:

m=-0.1;
h = 0.01;
f1 = @(x, y1, y2, y3) y2;
f2 = @(x, y1, y2, y3) y3;
f3 = @(x, y1, y2, y3) -0.5*(m+1)*y1*y3 +m*(y2^2)-m;
eta = 0:h:10;
x = 0:h:10;
y1(1) = 0;
y2(1) = 0;
y3(1) = 1;
for i = 1:(length(eta)-1)
a = h.*[f1(eta(i), y1(i), y2(i), y3(i)), f2(eta(i), y1(i), y2(i), y3(i)), f3(eta(i), y1(i), y2(i), y3(i))];
b = h.*[f1(eta(i), y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2), f2(eta(i)+h/2, y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2), f3(eta(i)+h/2, y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2)];
c = h.*[f1(eta(i), y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2), f2(eta(i)+h/2, y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2), f3(eta(i)+h/2, y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2)];
d = h.*[f1(eta(i), y1(i)+c(1), y2(i)+c(2), y3(i)+c(3)), f2(eta(i)+h, y1(i)+c(1), y2(i)+c(2), y3(i)+c(3)), f3(eta(i)+h, y1(i)+c(1), y2(i)+c(2), y3(i)+c(3))];
y3(i+1) = y3(i)+ 1/6*(a(3)+2*b(3)+2*c(3)+d(3));
y2(i+1) = y2(i)+ 1/6*(a(2)+2*b(2)+2*c(2)+d(2));
y1(i+1) = y1(i)+ 1/6*(a(1)+2*b(1)+2*c(1)+d(1));

i'm trying to solve it with finding the f''(0) initial value and than solve the equation, but it's seem not to work. could any one help me with that?

thank u

[Alan Stevens]

You have written

f3 = @(x, y1, y2, y3) -0.5*(m+1)*y1*y3 +m*(y2^2)-m;

but your mathematical equations suggest this should be:

f3 = @(x, y1, y2, y3) -0.5*(m+1)*y1*y2 +m*(y2^2)-m;

i.e. the first term in the RHS should have the y3 replaced by y2.

[eden meirovich]

are u sure? this this the function * the second derivative.

y1- the function f

y2 - the first derivative f'

y3 - the second derivative f ''

also, i'm not sure what do you ment with RHS

thank u!

[Alan Stevens]

My apologies. I had difficulty seeing the number of primes (') correctly! You are right it should be y1*y3 multiplying the 0.5(m+1) term. However, magnifying the whole equation I see that the term that is squared is also the second derivative also, so that shoud be y3^2 (unless my eyesight is really up the creek!).

Incidentally, by RHS I meant Right-Hand Side.

[eden meirovich]

hey! thanks for answring.

this is my fault now.this is the right equation:

m*f ' ^2 -0.5(m+1)*f*f'' = m+f'''

it is the first derivative

[John D'Errico]

You still have a problem. (Besides the fact you are writing your ODE solving code, instead of a better tool, like ODE45.) We see this statement:

my initial values are:

f(0) = f'(0) = 0

f(inf) =1

Yet, you have set initial values as:

y1(1) = 0;
y2(1) = 0;
y3(1) = 1;

So, are you trying to adjust the initial value for y''(0) as well as the value of m, such that the value y(inf) is 1? That is, do you want to use a shooting method here?

The problem is, there are TWO unknown pieces of information, thus y''(0) AND m.

[eden meirovich]

u can set m as u wish.

i trired using f ' ' (0) = 1 to find my next initial value f ' '(0) = ?

i look up on it online and find a way.

, do you want to use a shooting method here - i dont know the shooting method =\. so i can't write a code for that. i will be happy if u have one ready to use for me to try.

[John D'Errico]

Then read my answer. The solution exists, but has no limit as x-->inf.

[eden meirovich]

i have made a corection to my post.

f ' (inf) = 1

is that change staff with your answer?

i will also add that the solution to f ' is the one that i more care about.

[John D'Errico]

Sorry, but no. The solution to your problem appears to be oscillatory. It is no more true that we can say the limit of sin(x) exists as x-->inf than does the derivative have a limit.

You certainly will+ have no success with a numerical ODE solution, because as you can see from the plot I made, these ocillations in the solution seem to be of virtually constant magnitude. So your numerical solution will need to go out a really long way, if you would hope those oscillations will die down.

Then look at the analytical solution found by dsolve. Those oscillations will continue forever, because they are composed of terms like 9*cos(3*t))/80. A limit at infinity makes no sense here.

[eden meirovich]

hey. thank u for answring. im adding the equation from the book i'm using. u can see there how the solution look like. it's a numerical one so i don't quite sure how u find a analitc one.

[Alan Stevens]

I get reasonably sensible results using Matlab's own ode45 function, by guessing initial values and adjusting them to get f'(inf) = 1 (approx!) see below. One could presumably improve on this by using fzero, or bvp4c.

% m*f'^2 -0.5(m+1)*f*f'' = m+f'''
% f(0) = f'(0) = 0  f'(inf) = 1

m = 4;
% m = -0.09;
scale = (0.5*(m+1))^0.5;
L = 7/scale;
n = 1000;
xspan = 0:L/n:L;
% Now guess initial value of d2fdx2 needed to make dfdx(inf) = 0
d2fdx2_0 = 2.4057;  % for m = 4
% d2fdx2_0 = 10^-6;  % for m = -0.09
F0 = [0 0 d2fdx2_0];


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

dfdx = F(:,2);
xscale = scale*x;
plot(xscale, dfdx), grid
xlabel('(0.5(m+1))^{0.5}\eta'), ylabel('f''(\eta)')
legend(['m = ', num2str(m)])

function dfdx = rates(~,F,m)

         f = F(1);
         dfdx = F(2);
         d2fdx2 = F(3);
         d3fdx3 = m*dfdx^2 - 0.5*(m+1)*f*d2fdx2 - m;
        
         dfdx = [dfdx; d2fdx2; d3fdx3];
         
end

The resulting curves don't match those of your picture above exactly, but they are similar.

[eden meirovich]

thanks! i will give it a try. how to u guess f''(0)? just chaning numbers?

[Alan Stevens]

Yes!! (You could start by estimating the slope from your published graph).

[eden meirovich]

thanks! it actually works quite nice! hoped to find a good way to know the f''(0) value with out gussing. but i will work on if.

do u have tips how could i run a for loop and plot every m value in deifrent color and add a correct legend for that?

[Alan Stevens]

If you are guessing f''(0) "by hand" it will be difficult to automate doing several values of m at once. You probably need to get fzero involved to get an auromatic estimate of f''(0). Try using fzero for one value of m, then see about automating several at once.

[eden meirovich]

i have found servel values already:

m = [0 4 -0.09 1 1/3 1/2];

d2fdx2_0 = [0.33 2.4057 10^-6 1.2325 0.757 0.8995];

put them in a vector. just cant get it to plot it OK in a for loop with the ode45

[Alan Stevens]

Do something like

for i = 1:4
F0 = [0 0 d2fdx2_0(i)];
[x, F] = ode45(@rates, xspan, F0, [], m(i));
dfdx = F(:,2);
xscale = scale*x;
plot(xscale, dfdx), grid
hold on
end
xlabel('(0.5(m+1))^{0.5}\eta'), ylabel('f''(\eta)')

I've just made this up off the top of my head as I don't have time to test it right now! However, it should give you a "starter for ten".

Answer 1

[John D'Errico]

I think you have multiple problems here, now that you claim to have the correct equation. DSOLVE finds an analytical solution to your ODE.

syms y(x)
dy = diff(y,x,1);
ddy = diff(y,x,2);
dddy = diff(y,x,3);
m = -0.1;

ODE = m*dy(x)^2 -0.5*(m+1)*y(x)*ddy(x) == m + dddy(x)
ODE =
- diff(y(x), x)^2/10 - (9*y(x)*diff(y(x), x, x))/20 == diff(y(x), x, x, x) - 1/10

ysol = dsolve(ode,y(0) == 0,dy(0) == 0,ddy(0) == 1)
ysol =
(9*cos(3*t))/80 + (9*exp(-t))/20 + (7*sin(3*t))/80 - (5*cos(t))/16 + (9*sin(t))/16 - (cos(2*t)*(cos(3*t)/8 - sin(3*t)/8 + (3*cos(t))/8 + sin(t)/8))/2 - (tan(t/2)*cos(t)*(tan(t/2) - 3*tan(t/2)^2 - 6*tan(t/2)^3 + 3*tan(t/2)^4 + tan(t/2)^5 - tan(t/2)^6 + 1))/(tan(t/2)^2 + 1)^4
>> pretty(ysol)
                                                                     / cos(3 t)   sin(3 t)   3 cos(t)   sin(t) \
                                                            cos(2 t) | -------- - -------- + -------- + ------ |
cos(3 t) 9   9 exp(-t)   sin(3 t) 7   5 cos(t)   9 sin(t)            \     8          8          8         8   /
---------- + --------- + ---------- - -------- + -------- - ----------------------------------------------------
    80           20          80          16         16                                2

        / t \        /    / t \             / t \3        / t \4        / t \5      / t \6     \
     tan| - | cos(t) | tan| - | - 3 #1 - tan| - |  6 + tan| - |  3 + tan| - |  - tan| - |  + 1 |
        \ 2 /        \    \ 2 /             \ 2 /         \ 2 /         \ 2 /       \ 2 /      /
   - -------------------------------------------------------------------------------------------
                                                      4
                                              (#1 + 1)

where

            / t \2
   #1 == tan| - |
            \ 2 /
            

Now, I'm not very worried about the exact solution found, because you have no idea what y''(0) should be. Instead, I want you to look at the solution itself. The solution is a sum of trig terms. In fact, we can plot it.

fplot(ysol,[0,100])

Do you see the essential nature of the solution is an oscialltory function, that does not settle down? The limit as x approaches infinity will not exist. This is no different from asking what the limit of sin(x) is, as x approaches +inf. It has no limit.

Possibly the problem is still in the equation, that you have not told us the correct form, or that I got it wrong, or that m is the wrong value, or that y''(0) was chosen incorrectly. One thing you can see is that any value for y''(0) that I choose results in a similar infinitely oscillatory behavior, that is not converging to anything. Changing m also results in just a different oscillatory behavior. So my guess is you have given us the wrong equation for the ODE.

Answer 2

[Alan Stevens]

I think the following does all that you want. My implementation of the legend is cumbersome, but works!

% m*f'^2 -0.5(m+1)*f*f'' = m+f'''
% f(0) = f'(0) = 0  f'(inf) = 1

m = [0 4 -0.09 1 1/3 1/2];
d2fdx2_0 = [0.33 2.4057 10^-6 1.2325 0.757 0.8995];
% m = -0.09;
figure(1)
hold on
for i=1:length(m)
    scale = (0.5*(m(i)+1))^0.5;
    L = 7/scale;
    n = 1000;
    xspan = 0:L/n:L;
    % Now guess initial value of d2fdx2 needed to make dfdx(inf) = 0
    %d2fdx2_0 = 2.4057;  % for m = 4
    % d2fdx2_0 = 10^-6;  % for m = -0.09
    F0 = [0 0 d2fdx2_0(i)];

    [x, F] = ode45(@rates, xspan, F0, [], m(i));
    dfdx = F(:,2);
    xscale = scale*x;
    plot(xscale, dfdx)
      
end
grid
xlabel('(0.5(m+1))^{0.5}\eta'), ylabel('f''(\eta)')
lg1 = [' m = ', num2str(m(1))]; lg2 = [' m = ', num2str(m(2))];
lg3 = [' m = ', num2str(m(3))]; lg4 = [' m = ', num2str(m(4))];
lg5 = [' m = ', num2str(m(5))]; lg6 = [' m = ', num2str(m(6))];
legend(lg1,lg2,lg3,lg4,lg5,lg6);

function dfdx = rates(~,F,m)

         f = F(1);
         dfdx = F(2);
         d2fdx2 = F(3);
         d3fdx3 = m*dfdx^2 - 0.5*(m+1)*f*d2fdx2 - m;
        
         dfdx = [dfdx; d2fdx2; d3fdx3];
         
end

[eden meirovich]

thanks! i was able to do that, but i'm no trying to make every plit line in a diffrent color. do u know how i can do that?

[Alan Stevens]

Try "help plot" in the command window to see the various plot options.

[eden meirovich]

i know how to change a color of a singel line. just cant figure it how to change it doring a for loop.

[Alan Stevens]

Set up a color style list before the for loop, something like:

clr = ['k', 'c', 'm', 'r', 'b', 'g'];

then in the plot command:

plot(xscale, dfdx,clr(i))

I think that will work (though I haven't tried it!)

Answer 3

[MINATI PATRA]

@ eden meirovich

What about bvp 4c?