Почему мое уравнение (y '= (3 * x ^ 3-y) / (3 * x)) не смещено должным образом с моими функциями построения? - PullRequest
Новая
       50

Почему мое уравнение (y '= (3 * x ^ 3-y) / (3 * x)) не смещено должным образом с моими функциями построения?

0 голосов
/ 07 октября 2018

Я пытаюсь построить решение моего дифференциального уравнения, но я не могу получить правильный график с помощью этого метода.Согласно Desmos, моя функция должна выглядеть так:

my function

Это мой код: 

clear
syms Y(X)
ode = diff(Y,X) == (3.*X.^3-Y)./(3.*X);
cond = Y(1) == 5;
YSol(X) = dsolve(ode, cond)
[X,Y] = meshgrid(-5 : .2 : 5);
Z = @(X,Y)(3.*X.^3-Y)./(3.*X);
dirfield(Z,-10:0.5:10,-10:0.5:10)
hold on;
[Xs,Ys] = ode45(Z,[-1,10],5); plot(Xs,Ys)
hold off

Это функцияdirfield: 2 

function dirfield(f,tval,yval)
% dirfield(f, t1:dt:t2, y1:dy:y2)
%   
%   plot direction field for first order ODE y' = f(t,y)
%   using t-values from t1 to t2 with spacing of dt
%   using y-values from y1 to t2 with spacing of dy
%
%   f is an @ function, or an inline function,
%     or the name of an m-file with quotes.
%
% Example: y' = -y^2 + t
%   Show direction field for t in [-1,3], y in [-2,2], use
%   spacing of .2 for both t and y:
%
%   f = @(t,y) -y^2+t
%   dirfield(f, -1:.2:3, -2:.2:2)

[tm,ym]=meshgrid(tval,yval);
dt = tval(2) - tval(1); 
dy = yval(2) - yval(1);
fv = vectorize(f);
if isa(f,'function_handle')
  fv = eval(fv);
end
yp=feval(fv,tm,ym); 
s = 1./max(1/dt,abs(yp)./dy)*0.35;
h = ishold;
quiver(tval,yval,s,s.*yp,0,'.r'); hold on;
quiver(tval,yval,-s,-s.*yp,0,'.r');
if h
  hold on
else
  hold off
end
axis([tval(1)-dt/2,tval(end)+dt/2,yval(1)-dy/2,yval(end)+dy/2])

Кто-нибудь знает, что я сделал не так?Заранее спасибо.

РЕДАКТИРОВАТЬ:

График из Desmos имеет смешанные оси, но моей проблемы нет.Я не понимаю, почему я получаю это изображение:

image

вместо графика Desmos, повернутого против часовой стрелки (или что-то в этом роде).

Maybie, что мне нужно, это другая функция сюжета?

1 Ответ

0 голосов
/ 08 октября 2018

Это линейный ОДУ первой степени с особенностью при х = 0, точек сгиба от оси у нет.Действительно, решение получается путем интегрирования 

(x^(1/3)*y)' = x^(7/3)

так, чтобы 

(x^1/3)*y = 3/10*x^(10/3) + C  <==> y(x)  =  10/3*x^3 + C*x^(-1/3)

с некоторым выражением, содержащим большие числа для константы интегрирования C.Это также должно быть результатом, если вы напечатаете результат dsolve.

График, который вы включили, относится к другой проблеме.

Или график поменял оси, это график yx.Для начального условия y(1)=5 получается C=5/3 и график, который выглядит так, если зеркально отображать по диагонали.

Чтобы построить решение с обеих сторон от начального условия y(1)=5, вынужно интегрировать и построить дважды.Один раз для временного диапазона [1 5] и один раз для [1 0.1].Максимальная область решения равна (0, бесконечность), решения нет и, следовательно, нет графика для отрицательного x.Это было бы другое решение, например для начального условия y(-1)=-5.

...