NDSolve with arrays

svend

unread,

May 29, 2008, 8:44:15 AM5/29/08

to

hi all,
i want to solve a PDE by discretizing in one variable in an array and solve. the equations are of the following type:
n = 10;
x0 = 1;
x1 = 10;

play[x_] := IntegerPart[x]

eq1 := {D[f[j][x], x] == f[play[x]][x], f[j][x0] == j^2}

eqs = Flatten[Table[eq1, {j, 1, n}]]
fun = Flatten[Table[f[j], {j, 1, n}]];
tes = NDSolve[eqs, fun, {x, x0, x1}][[1]];

my problem is that the variable x, in which the equation is differential, also appears on the rhs in the "index" of my array. mathematica does not evaluate "play[x]" during NDSolve and just tells me that the rhs is not numerical. I tried "_NumericQ" and similar things but nothing seems. I would be very happy if someone could tell me how to resolve this.

thanks
svend

Oliver Ruebenkoenig

unread,

May 30, 2008, 2:57:17 AM5/30/08

to

How about,

n = 10;
x0 = 1;
x1 = 10;

play[x_] := IntegerPart[x]

eq1[j_] := {D[ToExpression["f" <> ToString[j]][x], x] ==
ToExpression["f" <> ToString[play[j]]][x],
ToExpression["f" <> ToString[j]][x0] == j^2}

eqs = Flatten[Table[eq1[j], {j, 1, n}]]
fun = Flatten[Table[ToExpression["f" <> ToString[j]], {j, 1, n}]]

tes = NDSolve[eqs, fun, {x, x0, x1}][[1]];

Oliver

Oliver Ruebenkoenig, <ruebenko AT uni-freiburg.de>

dh

unread,

Jun 1, 2008, 3:36:28 AM6/1/08

to

Hi Svend,

Name the 10 functions you are looking for: f1..f10. Note that all 10

functions have the same derivatives. Further, from t=m to t=m+1, the

derivative is fm[x]. We therefore solve the equations from t=1 to t=2

using f1 for the derivatives. Then from t=2 to t=3 using f2 etc.

Therefore we have an outermost Do loop that iterates over starting times.

To solve, we first setup a list with function names:funs =

Table[Symbol["f" <> ToString[j]], {j, 1, n}]

Then we define the equations, remembering that in step m the derivatives

are given by fm:

eq1[j_]:={D[funs[[j]][x],x]==funs[[x0]][x],funs[[j]][x0]==(funs[[j]][x0]/.tes[[x0]])};

eqs=Flatten[Table[eq1[j],{j,1,n}]];

finally we solve ad append the soulution to the solution list: tes:

AppendTo[tes, NDSolve[eqs, funs, {x, x0, x0 + 1}][[1]]];

Finally we need to assemple the different pieces of the 10 functions

into Piecewise functions:

Table[Piecewise@Table[{(funs[[i]]/.Transpose[Drop[tes,1]][[i,j]])[x],x<j+1},{j,1,n-1}]

,{i,1,10}];

Here is the whole code:

n=10;

funs=Table[Symbol["f"<>ToString[j]],{j,1,n}];

tes={Table[funs[[j]][1]->j^2,{j,1,n}]};

Do[

eq1[j_]:={D[funs[[j]][x],x]==funs[[x0]][x],funs[[j]][x0]==(funs[[j]][x0]/.tes[[x0]])};

eqs=Flatten[Table[eq1[j],{j,1,n}]];

AppendTo[tes,NDSolve[eqs,funs,{x,x0,x0+1}][[1]]];

,{x0,1,10}];

res=Table[Piecewise@Table[{(funs[[i]]/.Transpose[Drop[tes,1]][[i,j]])[x],x<j+1},{j,1,n-1}]

,{i,1,10}];

Plot[res,{x,1,4}]

hope this helps, Daniel

svend wrote:

> hi all,

> n = 10;

> x0 = 1;

> x1 = 10;

>

> thanks

> svend

>

--

Daniel Huber

Metrohm Ltd.

Oberdorfstr. 68

CH-9100 Herisau

Tel. +41 71 353 8585, Fax +41 71 353 8907

E-Mail:<mailto:d...@metrohm.com>

Internet:<http://www.metrohm.com>

Svend....@googlemail.com

unread,

Jun 1, 2008, 3:38:05 AM6/1/08

to

to be more precise,
play[x_, j_] :=If[ IntegerPart[x^2+j^2]>n, n, IntegerPart[x^2+j^2]>n],
i.e. the index on the rhs of my equation depends both on j and x.
any suggestions?
thanks
svend

On May 29, 2:44 pm, svend <dom...@tphys.uni-heidelberg.de> wrote:
> hi all,

> i want to solve a PDE by discretizing in one variable in an array and solv=

e. the equations are of the following type:
> n = 10;
> x0 = 1;
> x1 = 10;
>
> play[x_] := IntegerPart[x]
>
> eq1 := {D[f[j][x], x] == f[play[x]][x], f[j][x0] == j^2}
>
> eqs = Flatten[Table[eq1, {j, 1, n}]]
> fun = Flatten[Table[f[j], {j, 1, n}]];
> tes = NDSolve[eqs, fun, {x, x0, x1}][[1]];
>

> my problem is that the variable x, in which the equation is differential, =
also appears on the rhs in the "index" of my array. mathematica does not eva=
luate "play[x]" during NDSolve and just tells me that the rhs is not numeric=
al. I tried "_NumericQ" and similar things but nothing seems. I would be ver=

dh

unread,

Jun 2, 2008, 4:31:43 AM6/2/08

to

Hi Svend,

define the derivatives by:

deriv[j_Integer,x_?NumericQ,funvals_]:=funvals[[

If[IntegerPart[Sqrt[j^2+x^2]]>n,n,IntegerPart[Sqrt[j^2+x^2]]] ]];

and replace: D[funs[[j]][x],x]==funs[[x0]][x] by:

D[funs[[j]][x],x]== deriv[x0,x,Table[(funs[[i]])[x],{i,1,n}]]:

the whole code:

n=10;

funs=Table[Symbol["f"<>ToString[j]],{j,1,n}];

deriv[j_Integer,x_?NumericQ,funvals_]:=funvals[[

If[IntegerPart[Sqrt[j^2+x^2]]>n,n,IntegerPart[Sqrt[j^2+x^2]]] ]];

tes={Table[funs[[j]][1]->j^2,{j,1,n}]};

Do[eq1[j_]:={D[funs[[j]][x],x]==

deriv[x0,x,Table[(funs[[i]])[x],{i,1,n}]],funs[[j]][x0]==(funs[[j]][x0]/.tes[[x0]])};

eqs=Flatten[Table[eq1[j],{j,1,n}]];

AppendTo[tes,NDSolve[eqs,funs,{x,x0,x0+1}][[1]]];,{x0,1,10}];

res=Table[Piecewise@Table[{(funs[[i]]/.Transpose[Drop[tes,1]][[i,j]])[x],x<j+1},{j,1,n-1}],{i,1,10}];

Plot[res,{x,1,4}]

hope this helps, Daniel

Svend....@googlemail.com wrote:

> to be more precise,

> any suggestions?

> thanks

> svend

>

>> hi all,

>> n = 10;

>> x0 = 1;

>> x1 = 10;

>>

>> thanks

>> svend

>

--