# boundary problem

I have some conditions at t=0, like r1^2+r2^2=0.5 and x^2+y^2=6 and I want to find the right initial values in the past like t=-20, which reproduces this conditions. For example I have this system here.

alpha = 0.5
beta = 1

sol = NDSolve[{

x'[t] == -x[t]*a[t] - 3*x[t] + alpha*y[t]^2 +
1/2*beta*(r1[t]^2 - r2[t]^2),
y'[t] == -y[t]*a[t] - alpha*x[t]*y[t],
r1'[t] == -r1[t]*a[t] - 3/2*r1[t] + beta*x[t]*r1[t],
r2'[t] == -r2[t]*a[t] - 3/2*r2[t] - beta*x[t]*r2[t],
z'[t] == -z[t]*a[t] - 2*z[t],
x[-20] == 0, y[-20] == 10^(-12), r1[-20] == 0.008,
r2[-20] == 0.008, z[-20] == 0.999},
{x, y, r1, r2, z},
{t, 0, -20}];
a[t_] = -1/2*(3 - 3*y[t]^2 + 3*x[t]^2 + z[t]^2)


In c++ I would build something with a While loop. But in Mathematica I don't know and the documentation of Mathematica doesn't help me.

-
Does anyone know something about boundary value problems? – lambda1990 Jun 25 '13 at 12:28
I have no idea how to do it effectively. But you should define a[t_] before using it. – Silvia Jun 25 '13 at 15:00
There will almost certainly be multiple solutions corresponding to various combinations of {x[-20],y[-20],z[-20],r1[-20],r2[-20]}. If you have an idea of the range of these functions, or if you have other boundary conditions, the easiest way would be to include your criteria as conditions as well, i.e.:r1[0]^2+r2[0]^2==0.5 – Corey Kelly Jun 26 '13 at 16:12

First let us find all the solution as functions (I'm going to use ReplaceAll(\.) for that).

Clear[a,x, y, z, r1, r2]

alpha = 0.5;
beta = 1;

tmin = -20; tmax = 0;

a[t_] = -1/2*(3 - 3*y[t]^2 + 3*x[t]^2 + z[t]^2);

{{x, y, r1, r2, z}} = {x, y, r1, r2, z} /. NDSolve[{x'[t] == -x[t]*a[t] - 3*x[t] + alpha*y[t]^2 + 1/2*beta*(r1[t]^2 - r2[t]^2),
y'[t] == -y[t]*a[t] - alpha*x[t]*y[t],
r1'[t] == -r1[t]*a[t] - 3/2*r1[t] + beta*x[t]*r1[t],
r2'[t] == -r2[t]*a[t] - 3/2*r2[t] - beta*x[t]*r2[t],
z'[t] == -z[t]*a[t] - 2*z[t], x[-20] == 0, y[-20] == 10^(-12),
r1[-20] == 0.008, r2[-20] == 0.008, z[-20] == 0.999}, {x, y, r1,
r2, z}, {t, tmin,tmax}];


If you wish you can check your solutions by,

Plot[{x[t], y[t], z[t], r1[t], r2[t]}, {t,tmin,tmax}]


Now you want to make good guess for initial value and I think you want to do it iteratively (as you mention While). I am not sure what is convergence criterion, so I'll give you a general example. Say your initial values are m and n and your convergence condition is q = m+n < 5.65. Now you can use While (yes it is there in Mathematica as well)

m = 1; n = 1; (*initial guess*)
While[q = m + n; q < 5.65, Print[m, "  ", n, "  ", q]; m = m + 0.2; n = n + 0.3]


But I rather suggest you to use a Do loop.

Do[q = m + n; If[q > 5.5 && q < 5.65, Print[m, "  ", n, "  ", q]], {m, 1.3,4.1, .23}, {n, 1.7, 4.2, .21}]


Now you can scan the convergence within a particular range of initial values. Here I narrow down the condition just to show you that you can do it. You can use additional Break[] to stop the loop soon as you find your first set of initial values.

Say you convergence criterion is 0<x^2+y^2<0.001 at t=0. Then the structure will be like

Do[  Clear[a,x,y, ...
...
.
... x[-20]==xini, y[-20]==yini, ...

conv[t_]:=x[t]^2+y[t]^2;
If[conv[0] > 0.0 && conv[0] < 0.001, initval={xini,yini, ...}; Break[]];

,{xini,xini1,xini2,dxini}, {yini,yini1,yini2,dyini}, ...]


initval will record your answer. You can make a clever guess for a range of initial values and play with it whole day.

-