help solving differential equations

7 次查看（过去 30 天）
eden meirovich
eden meirovich 2020-8-1
回答： MINATI PATRA 2020-10-21
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 个评论
显示 15更早的评论
eden meirovich
eden meirovich 2020-8-2
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
Alan Stevens 2020-8-2
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".

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

John D'Errico
John D'Errico 2020-8-1
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 2020-8-3
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 个评论
显示 2更早的评论
eden meirovich
eden meirovich 2020-8-3
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
Alan Stevens 2020-8-3
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!)

请先登录，再进行评论。

MINATI PATRA
MINATI PATRA 2020-10-21
@ eden meirovich
What about bvp 4c?

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by