CODEWITHSUNDEEP
matlabdifferential-equationsnumerical-integrationode45
Ribbon Cable
Ribbon Cable

Reputation: 23

Matlab: Differential equation (ode45): Can I reverse tspan for better initial conditions?

I'm using ode45 to solve/plot a second-order differential equation in Matlab. My tspan is from 0 to 0.25. But the initial conditions near zero are ill-defined (slope goes to infinity, complex values). The conditions near 0.25 are well defined (both slope and value are zero).

Questions:

  1. Can I reverse tspan, and use the "final conditions" as initial conditions?

  2. Well, I know I can do it (see code below), and I get a plot that looks like what I expect, but is this a valid thing to do in general? Am I getting lucky in this one case?

  3. ode45 provides numerical solutions and is not exact. Am I likely to have larger error after reversing tspan?

Here's my code, which should run standalone:

function ReverseTspan()
% solve diff-eq backward from tspan end to tspan start using ode45()
%  - Good initial conditions at the end, but not start.
%  - Is this a valid thing to do?

% clean slate
clc; clear all; close all;

% tspan - reversed!
R = 0.25:-0.001:0;

% initial values
hinits=[0.0000001;0]; 

% solve
[R,v] = ode45(@equ7,R,hinits);

% strip imaginary values (can't plot 'em)
v(find(real(v)~=v)) = NaN;

% plot first column 
plot(R,v(:,1)); 


function vprime = equ7(R,v);
% Solve second order non-linear differential equation 7:
%   v''(R) + 2/R*v'(R) = K_plus/(R^2)*( v^(-1/2) - lamda_plus*(1-v)^(-1/2)
%
% Matlab ode45 only likes first derivatives, so let:  
%    v_1(R) = v(R)
%    v_2(R) = v'(R)
%
% And create a system of first order diff eqs:
%    v_1'(R) = v_2(R)
%    v_2'(R) = -2/R*v_2(R) + K_plus/(R^2)*( v_1(R)^(-1/2) - lamda_plus*(1-v_1(R))^(-1/2)
%
% Constant Parameters:
K_plus = 0.0859;
lambda_plus = 3.7;

% Build result in pieces for easier debugging of problematic terms
int1 = 1 - v(1);
int2 = int1^(-1/2);
int3 = v(1)^(-1/2);
int4 = K_plus/(R^2);
vprime2 = -2/R*v(2);
vprime2 = vprime2 + int4*( int3 - lambda_plus*(int2) );
vprime=[v(2); vprime2 ];

Upvotes: 2

Views: 989

Answers (1)

horchler
horchler

Reputation: 18484

Numerically, you're never going to be able to start at (or get to) R=0 with that ODE. NaN will be returned for the first step after your initial condition, thus corrupting any future results. It might be possible to rescale your system to get rid of the singularity, but that's a mathematics issue. You could also try starting your integration at a time just a bit larger than zero, provided you specify initial conditions that are valid (satisfy your ODE).

  1. Can I reverse tspan, and use the "final conditions" as initial conditions?

Yes. This is a standard procedure. In particular, it's often used to find unstable manifolds and unstable equilibria.

  1. Well, I know I can do it (see code below), and I get a plot that looks like what I expect, but is this a valid thing to do in general? Am I getting lucky in this one case?

It's hard to generalize across all possible systems. There are likely cases where this won't work for reasons of stability or numerical error.

  1. ode45 provides numerical solutions and is not exact. Am I likely to have larger error after reversing tspan?

I can't see why it would necessarily be larger, but again this will depend on the system and solver used. Try using other solvers (e.g., ode15s or ode113) and/or varying the absolute and relative tolerances via odeset.

Upvotes: 0

Related Questions