11 views (last 30 days)

Show older comments

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
on 2 Aug 2020

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".

John D'Errico
on 1 Aug 2020

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.

Alan Stevens
on 3 Aug 2020

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

Alan Stevens
on 3 Aug 2020

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!)

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

Start Hunting!