help solving differential equations

Question

eden meirovich 2020-8-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/573775-help-solving-differential-equations

回答： 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更早的评论隐藏 15更早的评论

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

请先登录，再进行评论。

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

Answer 1

John D'Errico 2020-8-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/573775-help-solving-differential-equations#answer_473839

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Alan Stevens 2020-8-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/573775-help-solving-differential-equations#answer_474361

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