我尝试使用 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 方式,使用\pstODESolve
(RKF45
方法):
\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}