MATLAB Answers

# help solving differential equations

4 views (last 30 days)
eden meirovich on 1 Aug 2020
Commented: Alan Stevens on 3 Aug 2020
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

#### 17 Comments

Alan Stevens on 2 Aug 2020
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
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
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".

Sign in to comment.

### Answers (2)

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.

#### 0 Comments

Sign in to comment.

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

#### 4 Comments

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

Sign in to comment.

### Community Treasure Hunt

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

Start Hunting!