help solving differential equations

11 views (last 30 days)
eden meirovich
eden meirovich on 1 Aug 2020
Answered: MINATI PATRA on 21 Oct 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
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 (3)

John D'Errico
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
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
Alan Stevens
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.


MINATI PATRA
MINATI PATRA on 21 Oct 2020
@ eden meirovich
What about bvp 4c?

Categories

Find more on Data Type Identification in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!