almost-autnomous differential equation

y' = y^2 − t. This differential equation has “stationary” solutions, but unlike with an autonomous equation, those stationary solutions are not horizontal. Vary the initial condition y(0) = c for a bit to try to get a sense of what the solutions look like. (Picking values between −3 and +3 should be good enough.)

Part A: There’s a value P such that, if y(0) < P, then the solution to the initial value problem decreases, while if y(0) > P, the solution to the initial value problem increases. Figure out what P is to two decimal places.

Attempted code:

    syms y(t);
    for c = -3:1:3
        fc = dsolve('Dy = y^2 - t' , y(0) == c);
    end

Not sure how I get it to print out an answer for p and get it to be for p to two decimal places.

Answer 1

[Star Strider]

You need to solve it symbolically, but then you can use the matlabFunction function to create an anonymous function from it:

syms y(t) c
fc(t,c) = dsolve(diff(y) == y^2 - t, y(0) == c);
fc_fcn = matlabFunction(fc)

fc_fcn = @(t,c) -(airy(3,t)+((sqrt(3.0).*gamma(2.0./3.0).^2.*3.0+3.0.^(2.0./3.0).*c.*pi.*2.0).*airy(1,t))./(gamma(2.0./3.0).^2.*3.0-3.0.^(1.0./6.0).*c.*pi.*2.0))./(airy(2,t)+((sqrt(3.0).*gamma(2.0./3.0).^2.*3.0+3.0.^(2.0./3.0).*c.*pi.*2.0).*airy(0,t))./(gamma(2.0./3.0).^2.*3.0-3.0.^(1.0./6.0).*c.*pi.*2.0));

You also need to decide on an appropriate range for ‘t’ if you haven’t already been given one. It might be easier to use a ribbon plot for this until you home in on the correct value for ‘c’.

[Bob]

Sorry I am still somewhat new to matlab. I tried putting your code into my code and that gave me multiple errors.

Error using sym/subsindex (line 732) Invalid indexing or function definition. When defining a function, ensure that the arguments are symbolic variables and the body of the function is a SYM expression. When indexing, the input must be numeric, logical, or ':'.

Error in sym/privsubsasgn (line 997) L_tilde2 = builtin('subsasgn',L_tilde,struct('type','()','subs',{varargin}),R_tilde);

Error in sym/subsasgn (line 834) C = privsubsasgn(L,R,inds{:});

Error in Project_2 (line 35) fc(t,c) = dsolve(diff(y) == y^2 - t, y(0) == c);

[Star Strider]

I’m not certain what you did when you put ‘[my] code into [your] code’. First, see if you can run my code alone as I posted it.

Otherwise, there could be version differences. I’m using R2016a, and the code I posted ran for me without error.

You can always just use the anonymous function form of the solution:

fc_fcn = @(t,c) -(airy(3,t)+((sqrt(3.0).*gamma(2.0./3.0).^2.*3.0+3.0.^(2.0./3.0).*c.*pi.*2.0).*airy(1,t))./(gamma(2.0./3.0).^2.*3.0-3.0.^(1.0./6.0).*c.*pi.*2.0))./(airy(2,t)+((sqrt(3.0).*gamma(2.0./3.0).^2.*3.0+3.0.^(2.0./3.0).*c.*pi.*2.0).*airy(0,t))./(gamma(2.0./3.0).^2.*3.0-3.0.^(1.0./6.0).*c.*pi.*2.0));

[Bob]

Ok now I am getting the same long answer as you for fc_fn but I do not think that is the answer I need. I am trying to find the value of p to two decimal places.

[Star Strider]

Apparently, ‘fc’ has some sort of interesting behaviour such that it neither increases nor decreases with time for some value of ‘c’. My guess is that you need to put ‘fc_fcn’ into a loop and find the value of ‘c’ that most closely produces the desired behaviour.

I’ll think about this to see if I can come up with a more efficient way to approach it, but ‘brute-forcing’ it may be the only way.

[Bob]

I saw an example that used:

for c = -4:2:4
    fc = dsolve('Dy + 2*t*y = 2', y(0) == c) ;
end

I dont know if that helps at all?

[Star Strider]

I don’t know what your function is supposed to do, only that you’re supposed to find a value of ‘c’ that does it.

This may give you some insight into the function’s behaviour:

t = linspace(0, 5, 50);
c = linspace(-3, 3, 50);
[T,C] = meshgrid(t,c);

fc_fcn = @(t,c) -(airy(3,t)+((sqrt(3.0).*gamma(2.0./3.0).^2.*3.0+3.0.^(2.0./3.0).*c.*pi.*2.0).*airy(1,t))./(gamma(2.0./3.0).^2.*3.0-3.0.^(1.0./6.0).*c.*pi.*2.0))./(airy(2,t)+((sqrt(3.0).*gamma(2.0./3.0).^2.*3.0+3.0.^(2.0./3.0).*c.*pi.*2.0).*airy(0,t))./(gamma(2.0./3.0).^2.*3.0-3.0.^(1.0./6.0).*c.*pi.*2.0));

figure(1)
meshc(T, C, fc_fcn(T,C))
grid on
axis([xlim    ylim    [-1  1]*1E+1])
xlabel('t')
ylabel('c')
zlabel('fc')

[Bob]

Does this provide me with a value of c though.

[Star Strider]

Not directly, but since I don’t understand what behaviour the particular value of ‘c’ is supposed to do, it should give you a way of determining the behaviour you’re looking for.

You can also plot the derivative (Jacobian) with meshc or contour if that would help:

t = linspace(0, 5, 50);
c = linspace(-3, 3, 50);
[T,C] = meshgrid(t,c);

fc_fcn = @(t,c) -(airy(3,t)+((sqrt(3.0).*gamma(2.0./3.0).^2.*3.0+3.0.^(2.0./3.0).*c.*pi.*2.0).*airy(1,t))./(gamma(2.0./3.0).^2.*3.0-3.0.^(1.0./6.0).*c.*pi.*2.0))./(airy(2,t)+((sqrt(3.0).*gamma(2.0./3.0).^2.*3.0+3.0.^(2.0./3.0).*c.*pi.*2.0).*airy(0,t))./(gamma(2.0./3.0).^2.*3.0-3.0.^(1.0./6.0).*c.*pi.*2.0));
dfc_fcn = gradient(fc_fcn(T, C));

figure(1)
contour(T, C, dfc_fcn, 50)
grid on
xlabel('t')
ylabel('c')