Как нарисовать трехмерное изображение: Plot3D NDSolve - PullRequest
3 голосов
/ 18 ноября 2011
m = 10; c = 2; k = 5; F = 12;

NDSolve[{m*x''[t] + c*x'[t] + (k*Sin[2*Pi*f*t])*x[t] == F*Sin[2*Pi*f*t], 
         x[0] == 0, x'[0] == 0}, x[t], {t, 0, 30}] 

{f, 0, 5} (0 =

Как нарисовать трехмерное изображение:

x = u (t, f)

............

Если f = 0,1,0,2, ... 5, Мы можем решить уравнение:

NDSolve[{m*x''[t] + c*x'[t] + (k*Sin[2*Pi*f*t])*x[t] == F*Sin[2*Pi*f*t], 
         x[0] == 0, x'[0] == 0}, x[t], {t, 0, 30}] 

x является функцией от t и f

...............

m = 10; c = 2; k = 5; F = 12;

f = 0.1

s = NDSolve[{m*x''[t] + c*x'[t] + (k*Sin[2*Pi*f*t])*x[t] == F*Sin[2*Pi*f*t], 
             x[0] == 0, x'[0] == 0}, x[t], {t, 0, 30}] 
Plot[Evaluate[x[t] /. s], {t, 0, 30}, PlotRange -> All]

f = 0,1 enter image description here

f = 0,2 enter image description here

f = 0,3 enter image description here

f = 5 enter image description here

Как нарисовать трехмерное изображение: x = u (t, f)

Ответы [ 3 ]

6 голосов
/ 18 ноября 2011

Здесь идет решение.

m = 10; c = 2; k = 5; F = 12;
NumberOfDiscrit$f = 20;(* Number of points you want to divide 0<=f<=5*)
NumberOfDiscrit$t = 100;(* Number of points you want to divide 0<=t<=30 *)
fValues = Range[0., 5., 5./(NumberOfDiscrit$f - 1)];
tValues = Range[0., 30., 30./(NumberOfDiscrit$t - 1)];
res = Map[(x /. 
  First@First@
    NDSolve[{m*x''[t] + c*x'[t] + (k*Sin[2*Pi*#*t])*x[t] == 
       F*Sin[2*Pi*#*t], x[0] == 0, x'[0] == 0}, x, {t, 0, 30}]) &,
fValues];
AllDat = Map[(#@tValues) &, res];
InterpolationDat = 
Flatten[Table[
Transpose@{tValues, 
  Table[fValues[[j]], {i, 1, NumberOfDiscrit$t}], 
  AllDat[[j]]}, {j, 1, NumberOfDiscrit$f}], 1];
Final3DFunction = Interpolation[InterpolationDat];
Plot3D[Final3DFunction[t, f], {t, 0, 30}, {f, 0, 5}, PlotRange -> All,
PlotPoints -> 60, MaxRecursion -> 3, Mesh -> None]

enter image description here

Вы можете использовать Manipulate для динамического изменения некоторых параметров. Кстати, вышеупомянутая трехмерная картинка может вводить в заблуждение, если принять f в качестве непрерывной переменной в u(t,f). Вы должны заметить, что численное решение кажется взорванным для асимптотических значений t>>30. Смотрите картинку ниже.

enter image description here enter image description here

Надеюсь, это поможет вам.

6 голосов
/ 18 ноября 2011

Вы также можете сделать что-то вроде этого

Clear[f]
m = 10; c = 2; k = 5; F = 12;

s = NDSolve[{m*Derivative[2, 0][x][t, f] + 
     c*Derivative[1, 0][x][t, f] + (k*Sin[2*Pi*f*t])*x[t, f] == F*Sin[2*Pi*f*t],
   x[0, f] == 0,
   Derivative[1, 0][x][0, f] == 0}, x, {t, 0, 30}, {f, 0, .2}]

Plot3D[Evaluate[x[t, f] /. s[[1]]], {t, 0, 30}, {f, 0, .2}, PlotRange -> All]

3d plot

3 голосов
/ 18 ноября 2011

Это должно сделать это.

m = 10; c = 2; k = 5; F = 12;


fun[f_?NumericQ] :=
 Module[
   {x, t}, 
   First[x /. 
     NDSolve[
      {m*x''[t] + c*x'[t] + (k*Sin[2*Pi*f*t])*x[t] == F*Sin[2*Pi*f*t],
       x[0] == 0, x'[0] == 0}, 
      x, {t, 0, 30}
     ]
   ]
 ]

ContourPlot[fun[f][t], {f, 0, 5}, {t, 0, 30}]

Важные моменты:

  • Шаблон _? NumericQ не позволяет fun быть оцененным для аргументов символов (например, fun[a]) и вызывать NDSolve::nlnum ошибок.

  • Поскольку NDSolve, по-видимому, не локализует свою переменную функции (t), нам нужно было сделать это вручную, используя Module, чтобы предотвратить конфликт между t, используемым в NDSolve, и один используется в ContourPlot. (Вы можете использовать переменную с другим именем в ContourPlot, но я думаю, что важно указать на это предостережение.)


Для значительного ускорения печати вы можете использовать памятка , как указал мистер Волшебник.

Clear[funMemo] (* very important!! *)

funMemo[f_?NumericQ] := 
 funMemo[f] = Module[{x, t}, 
   First[x /. 
     NDSolve[{m*x''[t] + c*x'[t] + (k*Sin[2*Pi*f*t])*x[t] == 
        F*Sin[2*Pi*f*t], x[0] == 0, x'[0] == 0}, x, {t, 0, 30}]]]

ContourPlot[funMemo[f][t], {f, 0, 5}, {t, 0, 30}] (* much faster than with fun *)

Если вы чувствуете себя предприимчивым и хотите более глубоко изучить Mathematica, вы можете улучшить это, ограничив объем памяти, который разрешено использовать в кэшированных определениях, , как я описал здесь .

Давайте определим вспомогательную функцию для включения запоминания:

SetAttributes[memo, HoldAll]
SetAttributes[memoStore, HoldFirst]
SetAttributes[memoVals, HoldFirst]

memoVals[_] = {};

memoStore[f_, x_] := 
 With[{vals = memoVals[f]}, 
  If[Length[vals] > 100, f /: memoStore[f, First[vals]] =.;
   memoVals[f] ^= Append[Rest[memoVals[f]], x], 
   memoVals[f] ^= Append[memoVals[f], x]];
  f /: memoStore[f, x] = f[x]]

memo[f_Symbol][x_?NumericQ] := memoStore[f, x]

Затем с использованием оригинальной незапамятной функции fun, вычерчивайте как

ContourPlot[memo[fun][f][t], {f, 0, 5}, {t, 0, 30}]
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...