# help solving differential equations

74 views (last 30 days)
eden meirovich on 1 Aug 2020 at 12:45
Commented: Alan Stevens on 3 Aug 2020 at 12:42
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 at 15:52
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 on 2 Aug 2020 at 16:05
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 on 2 Aug 2020 at 16:26
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 at 16:04
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 at 7:47
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

Show 1 older comment
Alan Stevens on 3 Aug 2020 at 11:29
Try "help plot" in the command window to see the various plot options.
eden meirovich on 3 Aug 2020 at 12:02
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 on 3 Aug 2020 at 12:42
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!)