Mathematica, численное решение системы дифуров.

[E@rthQu@ke]

Исходная задача: Имеется выражение для изменения вектора:
$\frac {d} {dz}\begin{bmatrix} S_1 \\ S_2 \\ S_3 \end{bmatrix}=\begin{bmatrix} S_1 \\ S_2 \\ S_3 \end{bmatrix}\times \begin{bmatrix} -1&0&0 \\ 0&0&0 \\ 0&0&0 \end{bmatrix}\begin{bmatrix} S_1 \\ S_2 \\ S_3 \end{bmatrix}+\begin{bmatrix} S_1 \\ S_2 \\ S_3 \end{bmatrix}\times \begin{bmatrix} 0&0&0 \\ 0&4&0 \\ 0&0&4 \end{bmatrix}\begin{bmatrix} P_1 \\ P_2 \\ P_3 \end{bmatrix}$
Проследить эволюцию системы на некоторых длинах z для разных состояний $S_1, S_2, S_3, P_1,P_2,P_3$, которые задаются, исходя из некоторых уравнений в цикле, либо могут быть подсчитаны предварительно и загнаны в массив, что несущественно, важно то, что считать нужно много раз. После этого обработать все получившиеся значения и вывести их как точку на сфере Пуанкаре($S_1, S_2,S_3$)
Я написал программу, которая мне численно считает(по Р-К) получившуюся систему уравнений, но я не уверен, насколько верно, потому что получаются не те данные, которые ожидаются.
И я не знаю, как организовать, чтобы начальные условия в численном расчёте менялись, а результаты куда-то записывались, чтобы их потом обработать и вывести. Разделить interpolating functions у меня не получается.
Вот код:

In[164]:= Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
ClearAll[s1,s2,s3,p1,p2,p3,S,P,l]
s0=1;
in1[theta_, fi_] = s0*Sin[theta]*Cos[fi];
in2[theta_, fi_] = s0*Sin[theta]*Sin[fi];
in3[theta_, fi_] = s0*Cos[theta];
p1=1;
p2=0;
p3=0;
S= {s1[z],s2[z],s3[z]}
P={p1[z],p2[z],p3[z]}
Jspm={{-js,0,0},{0,0,0},{0,0,0}};
Jxpm= {{0,0,0},{0,4*js,0},{0,0,4*js}};
rightpart=Simplify[(S\[Cross](Jspm.S)+S\[Cross](Jxpm.P))/js]
{4 0[z] (s2[z]-s3[z]),-s1[z] (4 0[z]+s3[z]),s1[z] (4 0[z]+s2[z])}
sol=NDSolve[{s1'[l]==4 p3 s2[l]-4 p2 s3[l], s2'[l]==-s1[l] (4 p3+s3[l]),s3'[l]==s1[l] (4 p2+s2[l]),s1[0]==0.19, s2[0]==0.19 ,s3[0]==0.96},{s1,s2,s3},{l,0,20}]
{{s1->InterpolatingFunction[{{0.,20.}},<>],s2->InterpolatingFunction[{{0.,20.}},<>],s3->InterpolatingFunction[{{0.,20.}},<>]}}
{Red,Green,Blue}]
out1=sol[[1]][[1]][[2]]
InterpolatingFunction[{{0.,20.}},<>]
out2=sol[[1]][[2]][[2]]
InterpolatingFunction[{{0.,20.}},<>]
out3=sol[[1]][[3]][[2]]
InterpolatingFunction[{{0.,20.}},<>]
Plot[Evaluate[out3],{l,0,20}]

Прошу совета опытных пользователей, заранее спасибо.

[Vince Diesel]

Код лучше копировать сюда без In и Out, чтобы его можно было скопировать и исполнить. Начальные условия можно менять в цикле. А результаты собирать в массив, можно с помощью AppendTo (хотя для больших массивов это не лучший способ):

h1[x_]=...
...
results={};
For[k=1,k<100,k++,
sol=Dsolve[...,s1[0]=h1[k]}...];
AppendTo[results,{s1,s2,s3}/.sol]
]

[E@rthQu@ke]

Исправил, кажется, спасибо. У меня не так и много получится, около 100 точек. За подсказку спасибо.

[Vince Diesel]

Если названия функций к нулю приравнять,

p1=1;
p2=0;
p3=0;

дальше чушь пойдет.
Может, имелось в виду

p1[z_]==1;
p2[z_]=0;
p3[z_]=0;

[E@rthQu@ke]

Vince Diesel в сообщении #579348 писал(а):

Если названия функций к нулю приравнять,

p1=1;
p2=0;
p3=0;

дальше чушь пойдет.
Может, имелось в виду

p1[z_]==1;
p2[z_]=0;
p3[z_]=0;

Значения p1, p2, p3 задаются вручную и пока не должны быть функциями. Просто вместо p2 должно подставляться 0 и т.д. При такой формулировке написанное верно?
Подскажите ещё пожалуйста, есть ли функция постройки точек на сфере, или мне строить сферу в Plot, а там уже множество точек, записанных в декартовых координатах? Либо перейти в сферическую ск и построить уже там, используя встроенную функцию SphericalPlot3d?

[Vince Diesel]

Если у p не предполагается аргументов, зачем писать далее

P={p1[z],p2[z],p3[z]}

Есть примитив Sphere. Для массива точек в декартовых координатах можно совмещать примитивы командой Graphics3D:

points = {{0, 0, 1}, {0, 1, 0}, {0, 0, 1}};
Graphics3D[{Sphere[{0, 0, 0}, 1], {Black, Point /@ points}}]

[E@rthQu@ke]

Vince Diesel

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
ClearAll[s0,s1,s2,s3,p1,p2,p3,S,P,l,sol,sum,results]
s0=1;
NP=100;
n=1;
in1[k_, n_] := s0*Sin[(k-1)*Pi/(NP-1)]*Cos[2(n-1)*Pi/(NP-1)]
in2[k_, n_] := s0*Sin[(k-1)*Pi/(NP-1)]*Sin[2(n-1)*Pi/(NP-1)]
in3[k_,n_] := s0*Cos[(k-1)*Pi/(NP-1)]
results={};
S= {s1[z],s2[z],s3[z]};
P={p1[z],p2[z],p3[z]};
For[k=1,k<101,k++,sol=NDSolve[
{s1'[l]==4 p3[l] s2[l]-4 p2[l] s3[l], s2'[l]==-s1[l] (4 p3[l]+s3[l]),s3'[l]==s1[l] (4 p2[l]+s2[l]),p1'[l]==4 s3[l] p2[l]-4 s2[l] p3[l], p2'[l]==-p1[l] (4 s3[l]+p3[l]),p3'[l]==p1[l] (4 s2[l]+p2[l]),
s1[0]==in1[k,k], 
s2[0]==in2[k,n],
s3[0]==in3[k,n],
p1[0]==1,
p2[0]==0,
p3[0]==0},{s1,s2,s3,p1,p2,p3},{l,0,5}];AppendTo[results,{s1,s2,s3,p1,p2,p3}/.sol]]
results[[1]]
Out[799]= {{InterpolatingFunction[{{0.,5.}},<>],InterpolatingFunction[{{0.,5.}},<>],InterpolatingFunction[{{0.,5.}},<>],InterpolatingFunction[{{0.,5.}},<>],InterpolatingFunction[{{0.,5.}},<>],InterpolatingFunction[{{0.,5.}},<>]}}
sol
Out[800]= {{s1->InterpolatingFunction[{{0.,5.}},<>],s2->InterpolatingFunction[{{0.,5.}},<>],s3->InterpolatingFunction[{{0.,5.}},<>],p1->InterpolatingFunction[{{0.,5.}},<>],p2->InterpolatingFunction[{{0.,5.}},<>],p3->InterpolatingFunction[{{0.,5.}},<>]}}

Подскажите пожалуйста, как вытащить теперь из этого массива значения функций в определённых точках? Просто Plot[Evaluate] работает с результатом sol, но не работает с частью выходного result.
В конце мне потребуется посчитать среднее по каждому элементу, например в какой-то точке. Также, с sol подстановка работает, а с results нет.

[E@rthQu@ke]

Вот для sol можно выполнить

In[81]:= s1[5] /. %

Out[81]= {0.84971}

In[82]:= s2[3] /. %%

Out[82]= {-0.724476}

, а для results нельзя. Как в results добавить эти самые стрелочки?) Либо вытаскивать как-то иначе?

P.S. Подумал и сотворил

q1 = results[[1]][[1]][[1]]
In[267]:= q1[5]
Out[267]= 0.84971

[Vince Diesel]

Можно так:

f[t_] = results /. InterpolatingFunction[a$___] -> InterpolatingFunction[a$][t]
Plot[f[t][[1;;3]], {t, 0, 1}]