在 Asymptote 中求解 ODE

在 Asymptote 中求解 ODE

我尝试使用 Asymptote 根据曲率和扭转绘制空间曲线,但在这种情况下无法使 ODE 求解器工作:

import ode;

real kappa(real t){return 1;}
real tau(real t){return 0.1 ;}

real[] func(real t,real[] x){
  real[] y=x;
  real a=kappa(t);
  real b=tau(t);
  y[0:2]=x[3:5];
  y[3:5] = a*x[6:8];
  y[6:8] = -a*x[3:5]+b*x[9:11];
  y[9:11]=-b*x[6:8];
  return y; 
}

real[] x0={0,0,0,1,0,0,0,1,0,0,0,1};
Solution S=integrate(x0,func,0,20,h=0.1,Euler);
S.y=transpose(S.y);

write(S.y[0]);

现在S.y[0]应该包含空间曲线的 x 坐标,但它只返回零。我测试了各种选项,但没有变化。

我已经在另一个程序(Gauss)中测试了该代码,并且它在那里可以正常工作。所以问题一定是我的 Asymptote 错误。我做错了什么?

答案1

将函数的第一行从 改为real[] y = x;real[] y = copy(x);解释:渐近线数组的工作方式与 Java 数组类似。具体来说,说real[] y = x被翻译为“让y成为数组 的另一个名称x”。然后你的更改操作y实际上是修改x,这显然把事情搞砸了。

import ode;

real kappa(real t){return 1;}
real tau(real t){return 0.1 ;}

real[] func(real t,real[] x){
  real[] y=copy(x);
  real a=kappa(t);
  real b=tau(t);
  y[0:2]=x[3:5];
  y[3:5] = a*x[6:8];
  y[6:8] = -a*x[3:5]+b*x[9:11];
  y[9:11]=-b*x[6:8];
  return y; 
}

real[] x0={0,0,0,1,0,0,0,1,0,0,0,1};
Solution S=integrate(x0,func,0,20,h=0.1,Euler);
S.y=transpose(S.y);

write(S.y[0]);

输出:

0:  0
1:  0.1
2:  0.2
3:  0.299
4:  0.396
5:  0.4900101
6:  0.5800606
7:  0.66521199799
8:  0.74456478392
9:  0.817268928670301
10: 0.88253296910301
11: 0.939632593356149
12: 0.987918634621347
13: 1.02682438669365
14: 1.05587216055627
15: 1.07467900811329
16: 1.0829615468487
17: 1.08053982760215
18: 1.06734019673244
19: 1.04339711360394
20: 1.00885389448845
21: 0.963962364525555
22: 0.909081410228328
23: 0.844674436049393
24: 0.771305739627151
25: 0.689635831400811
26: 0.600415735204237
27: 0.504480317110514
28: 0.402740700091229
29: 0.296175831869127
30: 0.185823282576104
31: 0.0727693573812032
32: -0.0418613829677167
33: -0.156917093826187
34: -0.271230004716683
35: -0.383628052959534
36: -0.492946678154748
37: -0.598040660015069
38: -0.697795880426028
39: -0.791140890170835
40: -0.877058161523339
41: -0.954594909885117
42: -1.02287337081551
43: -1.08110042315606
44: -1.12857645445138
45: -1.16470337147282
46: -1.1889916663043
47: -1.20106645708391
48: -1.20067243203384
49: -1.18767763576722
50: -1.16207604793707
51: -1.12398891598566
52: -1.07366481595009
53: -1.01147842786306
54: -0.937928025134943
55: -0.853631690285405
56: -0.759322282382004
57: -0.65584119440672
58: -0.544130951379378
59: -0.425226712288528
60: -0.300246750588747
61: -0.170381999094851
62: -0.0368847554200089
63: 0.0989433464456911
64: 0.235763984341133
65: 0.372215294437474
66: 0.506925388291969
67: 0.638526107672646
68: 0.765666880631573
69: 0.887028539903007
70: 1.00133696368006
71: 1.1073763992041
72: 1.20400233139496
73: 1.29015376195387
74: 1.36486476896568
75: 1.42727522298176
76: 1.47664054283129
77: 1.51234038292871
78: 1.53388615354352
79: 1.54092728629076
80: 1.5332561688872
81: 1.51081168589211
82: 1.47368131559126
83: 1.4221017472629
84: 1.35645799764707
85: 1.27728102038388
86: 1.18524381734446
87: 1.08115607599916
88: 0.965957372098677
89: 0.840708991830605
90: 0.706584442104337
91: 0.56485873156058
92: 0.416896518151569
93: 0.264139231553797
94: 0.108091290122693
95: -0.0496944575471042
96: -0.20763192724714
97: -0.364117482925951
98: -0.517545956139565
99: -0.666326842775628
100:    -0.80890051525468
101:    -0.943754286621699
102:    -1.06943816278465
103:    -1.18458012065271
104:    -1.28790075307666
105:    -1.37822712628201
106:    -1.45450570188128
107:    -1.51581418350511
108:    -1.56137215753994
109:    -1.59055040832136
110:    -1.60287880031163
111:    -1.59805263317786
112:    -1.57593739016094
113:    -1.53657181554892
114:    -1.48016927329628
115:    -1.40711735570659
116:    -1.31797572845661
117:    -1.21347221591399
118:    -1.09449714851397
119:    -0.962096011733206
120:    -0.817460453752456
121:    -0.6619177260532
122:    -0.496918647771045
123:    -0.324024200455752
124:    -0.144890874797971
125:    0.0387450952844123
126:    0.225084463202255
127:    0.412282505657726
128:    0.598467195034853
129:    0.781757831104838
130:    0.96028394850497
131:    1.13220431181094
132:    1.29572580723702
133:    1.4491220391138
134:    1.59075144033749
135:    1.71907470896613
136:    1.83267138804736
137:    1.93025541256803
138:    2.01068945606943
139:    2.07299791990389
140:    2.11637842023204
141:    2.14021164156917
142:    2.14406944086195
143:    2.12772110257489
144:    2.09113766293512
145:    2.03449424015934
146:    1.95817032698792
147:    1.86274802199089
148:    1.74900819669128
149:    1.61792461636956
150:    1.47065605326126
151:    1.30853645152763
152:    1.13306322365606
153:    0.945883777624059
154:    0.748780393033133
155:    0.543653582288204
156:    0.332504089573641
157:    0.117413695677966
158:    -0.0994749895224016
159:    -0.315979553049117
160:    -0.529899419181656
161:    -0.7390378918284
162:    -0.941224380341408
163:    -1.13433658614695
164:    -1.31632242571104
165:    -1.48522146575505
166:    -1.63918564929938
167:    -1.77649909603958
168:    -1.89559676772186
169:    -1.99508179853414
170:    -2.07374130199243
171:    -2.13056047928552
172:    -2.16473486942849
173:    -2.17568059873067
174:    -2.16304250585163
175:    -2.12670003892541
176:    -2.06677084269009
177:    -1.98361197606162
178:    -1.87781872392198
179:    -1.75022099082411
180:    -1.60187728861464
181:    -1.43406635439784
182:    -1.24827645956604
183:    -1.04619249455481
184:    -0.829680937301971
185:    -0.600772835854127
186:    -0.361644956939533
187:    -0.114599272382812
188:    0.137959026238998
189:    0.393544777511875
190:    0.649617142619737
191:    0.90360470547473
192:    1.15293113518926
193:    1.3950411573785
194:    1.62742657510233
195:    1.84765207713663
196:    2.0533805707624
197:    2.24239777840909
198:    2.41263584229108
199:    2.56219568861114
200:    2.68936791292405

答案2

PSTricks 方式,使用\pstODESolveRKF45方法):

\documentclass{article}

\usepackage{pst-ode}

\pstVerb{
  /kappa {pop 1} def
  /tau {pop 0.1} def
  /coefA {t kappa} def
  /coefB {t tau} def
}

\def\odeRHS{
  x[3] |
  x[4] |
  x[5] |
  coefA*x[6] |
  coefA*x[7] |
  coefA*x[8] |
  -coefA*x[3] + coefB*x[9] |
  -coefA*x[4] + coefB*x[10] |
  -coefA*x[5] + coefB*x[11] |
  -coefB*x[6] |
  -coefB*x[7] |
  -coefB*x[8]
}

\def\init{ 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 }

\pstVerb{
  /printArr {
    0 exch
    {
      exch dup 20 string cvs print 1 add
      (:\string\t) print
      exch 256 string cvs print
      (\string\n) print
    } forall pop
  } def
}

%RKF45
\pstODEsolve[algebraicAll]{SOL}{x[0]}{0}{20}{201}{\init}{\odeRHS}

\pstVerb{ [ SOL ] printArr } %write result to terminal

\begin{document}
dummy text
\end{document}

相关内容