Ordinary Differential Equations part 2

In solving ordinary differential equations, variations of problem from the general ode problem that we introduced in part 1 may appear. Here, we will show two such variations :
1. Consider the first example that we introduced in part 1. If the problem is changed to be as follows :

Solve the equation for t from 0 to 5 :

dy/dx = -2x + y

Initial condition : xt=2.5 = 1

We can see that the time of the initial condition is t = 2.5 but not t = 0 but the range of the solution we need to find is still setted from 0 to 5. So, how can we solve this problem ?
In this case, we need to separate the problem into two parts. The first part is from t=2.5 to t=5. The second part is from t=2.5 to t=0. Here is how we do it :

Input 1 : r=rkode('-2*t0+x0','1',2.5,5,0.5)
Ans. set r = 2.5, 3, 3.5, 4, 4.5, 5, 1, -1.891, -7.304, -16.876, -33.304, -61.032

Input 2 : d=r(1)
 

For the second part :

Input 3 : s=rkode('-2*t0+x0','1',2.5,0,-0.5)
Ans. set r = 2.5, 2, 1.5, 1, 0.5, 0, 1, 2.359, 2.791, 2.66, 2.187, 1.507

See the the interval value is now a negative value -0.5.

Input 4 : f=s(1)
 

We can then reverse f(0) using the function 'rev'and join it with d(0) to get the solution of the whole range.

Input 5 : c=rev(f(0))
 

Input 6 : h=join(c,d(0))
Ans. set h = 1.507, 2.187, 2.66, 2.791, 2.359, 1, 1, -1.891, -7.304, -16.876, -33.304, -61.032

To get a plot of the solution, we also needed to reverse and combine the t values :

Input 7 : c=rev(s(0))
 

Input 8 : k=join(c,r(0))
Ans. set k = 0, 0.5, 1, 1.5, 2, 2.5, 2.5, 3, 3.5, 4, 4.5, 5

Input 9 : plot('line','k:h',0,5,0.5)
 

The plot is shown below :

The exact solution of this problem is the formula below :
y = 2 - 6e^(x-2.5) + 2x
We can check our numerical solution h by using k as x and calculate y :

Input 4 : 2-6*exp(k-2.5)+2*k
Ans. 1.507, 2.188, 2.661, 2.793, 2.361, 1, 1, -1.892, -7.31, -16.89, -33.334, -61.095

2. The second variation from the general ode problem is also about the change in initial condition. Consider the following example

Solve the equation for x from 0 to 10 :

d2y/dx2 = √(1 + (dy/dx)2)

y(0) = 2, y(10) = 0

For this problem, only the values of the boundary are given but the initial value of dy/dx is unknown. These kind of problem is called 'boundary value problem'. So, how can we solve this problem ?
In this case, we can try to use a method called 'shooting'. The idea is quite simple. We try different initial values of dy/dx and see if the resulted value of y(10) is 0. If it is, then we had done.
First, we rename variables and reduce the problem into two first order odes :

Solve the equation for t0 from 0 to 10 :

dx0/dt0 = x1
dx1/dt0 = √(1 + (x1)2)

x0t0=0 = 2, x0t0=10 = 0

We set x as the initial value of dx0/dt0_t0=0 and try to find x0(10). Below are the inputs that can do the search :

Input 10 : a='x1|sqrt(1+x1*x1)'
 

Input 11 : b='2|x'
 

Input 12 : c='r=euode(a,b,0,10,0.2)|d=r(1)|f=d(0)'
 

Here, we use the function 'euode' instead of 'rkode' since 'euode' is a much faster method.

Input 13 : h='f(size(f)-1)'
 

With the above inputs, we can now do the shooting. For example, we set x=0 to see what is y(10) :

Input 14 : table(c,0,x,h:)
Ans. 4848.1

y(10) is found to be 4848.1.
Now, we try x = 50,

Input 15 : table(c,50,x,h:)
Ans. 455023.5

For x = -50,

Input 16 : table(c,-50,x,h:)
Ans. 57.33

For x = -80,

Input 17 : table(c,-80,x,h:)
Ans. -6.35

So, it seems that zero should exists between -50 to -80. We can continue to try to input different x to find the output. Or, we can actually do a zero search between -50 to -80 as follows with error level set to 0.1 (default is 0.001) :

Input 18 : zeros('table(c,s,x,h:)',-80,-50,'s',0.1)
Ans. -76.4375

It takes some time for this search step and finally we get the answer -76.4375. When we put the value x= -76.4375, y(10) would be -0.061.

Input 19 : table(c,-76.4375,x,h:)
Ans. -0.061

It can be seen that when the initial dy/dx is set to -76.4375 at x=0, then we get y(10) = -0.061 which is quite close to 0 and so we treat the calculated r and d from c as our solution.
To get a plot of the solution :

Input 20 : plot('line','r(0):d(0)',0,10,0.2)
 

The resulting graph is :