0
$\begingroup$

I want to solve the following Cauchy problem \begin{equation} y' = y^2 + \frac{t^4 - 6t^3 + 12t^2 - 14t + 9}{(1+t)^2} \end{equation}

with initial condition: $y(0) = 2$ for $t \in [0,1.6]$ using a 3 steps Adam-Moulton method implemented in Matlab as follows:

function [t,u] = AM3(fun,t0,T,y0,N)
  h = (T-t0)/N; % integration step
  t = t0:h:T; %time mesh
  [t(1:3),u(1:3)] = RK4(fun,t0,t0+2*h,y0,2); % first 3 steps

  for n=1:N-2
    fn = feval(fun,t(n),u(n));
    fn1 = feval(fun,t(n+1),u(n+1));
    fn2 = feval(fun,t(n+2),u(n+2));
    F = @(z) z - u(n+2) - (h/24)*(9*feval(fun,t(n+3),z) + 19*fn2 - 5*fn1 + fn);
    z0 = u(n+1);
    z1 = u(n+2);
    dF = @(z) 1 - (h/24)*9*(2*z);
    u(n+3) = newton(F,dF,z0);
  end
endfunction

N is the number of integration nodes used (i.e. the size of the time mesh). The function returns a discretized solution of the problem in the array $u$. Since the method is implicit at every iteration I use Newton's method written as

function y = newton(f,df,x)
 fx = f(x);
 tol = 1e-10;
 itermax = 1e3;
 iter = 0;
 while abs(fx) > tol && iter < itermax
    x = x - f(x)/df(x);
    fx = f(x);
    iter++;
 endwhile
 y=x;

My problem is that it seems that the equation is very ill-conditioned since even with a tiny integration step $h$ what I obtain is something like this: enter image description here the problem persisting with much bigger $N$.

Despite the equation being pretty ugly it has en exact solution: $y(t) = \frac{(t-1)(t-2)}{t+1}$ which is relatively simple.

Can you help me please getting around this? I double checked the algorithms and they seems fine; I also tried to change Newton's method with some other (secant, fixed point, etc) but no novelty occured.

I also tried changing altoghether the method used with an Adam-Bashforth or Runge-Kutta one but to no avail.

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ Did you try to run it for the simple test problem $y'= -5y, y(0)=1$ ? Does it have the right order for that equation? If not, your implementation is wrong $\endgroup$
    – VoB
    Jan 7, 2021 at 21:11
  • 2
    $\begingroup$ If you claim that $y(t) = \frac{(t-1)(t-2)}{t+1}$ has to be a solution, then the equation should have a minus, not a plus: $$\frac{dy}{dt} \, = \, y^2 \, - \, \frac{t^4 - 6t^3 +12t^2 - 14t + 9}{(t+1)^2}$$ $\endgroup$
    – Futurologist
    Jan 8, 2021 at 5:10

1 Answer 1

Reset to default
2
$\begingroup$

You can de-singularize Riccati equations $y'=y^2+a$ by setting $y=-\frac{u'}{u}$ to get $u''+au=0$. If $a$ is continuous on the integration interval, the solution will be well-behaved on that interval. Solving with $u(0)=1$ and $u'(0)=-y(0)=-2$ gives the plot

enter image description here

with a root of $u$ and thus a singularity of $y$ at about $t=0.35230$. In any numerical integration, all that follows after that is to be considered random noise.

In the task as you obviously intended it, you have an equation $y'=y^2+p'-p^2$, that is, $a=p'-p^2$, for some given solution function $p$. The 4th degree term in the nominator comes from the second term, so needs to have a negative coefficient. With this corrected sign in the equation the solution for the transformed equation looks like

enter image description here

so the denominator $u$ is a solid distance away from zero. Solving the equation directly also confirms this. There is no divergence even under heavy distortion of the initial value

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.