# Numerical solution of ill-conditioned differential equation

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

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

• 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 – VoB Jan 7 at 21:11
• 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}$$ – Futurologist Jan 8 at 5:10

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