intge1te.pas 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498
  1. program intge1te;
  2. uses
  3. typ,
  4. spe,
  5. int;
  6. const
  7. e = 2.71828182845905;
  8. fnames = 'KI A0 A1 A2 A3 A4 SS SL SE V1 V2 ';
  9. ogs: array[1..11] of ArbFloat = (0, 0, 1, e, 0, 1, 0, 1, 1, 1, 1);
  10. integraaltekst: array[1..11, 1..5] of string[60] =
  11. ((' ì ',
  12. ' ô -àcosh(x) ',
  13. {k0} ' ³ e dx = k0(à), mits à > 0. ',
  14. ' õ ',
  15. '0 '),
  16. (' ì ',
  17. ' ô sin x àcos x ',
  18. {a0} ' ³ ------- + ---------- dx = 1, mits à>0 ',
  19. ' õ (x+1)à (x+1)(à+1) ',
  20. '0 '),
  21. (' ì ',
  22. ' ô 1 ',
  23. {a1} ' ³ ---- dx = 1/(à-1), mits à>1 ',
  24. ' õ xà ',
  25. '1 '),
  26. (' ì ',
  27. ' ô dx ',
  28. {a2} ' ³ --------- = 1/(à-1), mits à>1 ',
  29. ' õ x.ln(x)à ',
  30. 'e '),
  31. (' ì ',
  32. ' ô Ú àùxàùsin(xà) cos(xà)¿ ',
  33. {a3} ' ³ ³ -------------- + ---------³ dx = 1, mits à>0 ',
  34. ' õ À x(x+1) (x+1)ý Ù ',
  35. '0 '),
  36. (' ì ',
  37. ' ô Ú 2sin(«ãùxà) xàùcos(«ãùxà) ¿ ',
  38. {a4} ' ³ ³-------------- + ãà-----------------³ dx = 1, mits àò0',
  39. ' õ À (x+1)ý x(x+1) Ù ',
  40. '1 '),
  41. (' ss(n)=2*(n+1)(à-1)/n (n=1,2,3...), àò0 ',
  42. {ss} ' ss(x)=0 als min(|n-x|) ò 0.5/(n+1)à ',
  43. ' ss lineair interpoleren als min(|n-x|) ó 0.5/(n+1)à ',
  44. ' int. 0:ì = ä [1:ì] 1/(n(n+1)) = 1 ',
  45. ' '),
  46. (' ì ',
  47. ' ô sin(ln(x)) 1 ',
  48. {sl} ' ³ --------- dx = ---------, mits à>1 ',
  49. ' õ xà (à-1)ý+1 ',
  50. '1 '),
  51. (' ì ',
  52. ' ô sin(xà)-à.x(à-1).cos(xà) sin(1) ',
  53. {se} ' ³ --------------------------- dx = ------ ',
  54. ' õ ex e ',
  55. '1 '),
  56. (' ì ',
  57. ' ô à.|x|(à-1) ',
  58. {v1} ' ³ ---------------- dx = 1, mits à > 0 ',
  59. ' õ ã.(|x|(2à) + 1) ',
  60. '-ì '),
  61. (' ì 0 ì ',
  62. ' ô ô àx ô -x/à ',
  63. {v2} ' ³ v2(x)dx = ³ e dx + ³ e dx = à + 1/à, mits à > 0',
  64. ' õ õ õ ',
  65. '-ì -ì 0 '));
  66. var
  67. alfa, ond, inte, int1: ArbFloat;
  68. u, i: ArbInt;
  69. s: string;
  70. q: char;
  71. f: rfunc1r;
  72. scale: boolean;
  73. function Ki(x: ArbFloat): ArbFloat;
  74. var
  75. kk: ArbFloat;
  76. begin
  77. if abs(x) < ln(100 / alfa) then
  78. kk := Exp(-alfa * Specoh(x))
  79. else
  80. kk := 0;
  81. if scale then
  82. ki := kk / int1
  83. else
  84. ki := kk;
  85. end;
  86. function uki(u: ArbFloat): ArbFloat; {u=1/(x+1), of x=1/u-1}
  87. begin
  88. if u > 0 then
  89. uki := ki((1 - u) / u) / sqr(u)
  90. else
  91. uki := 0;
  92. end;
  93. function a0(x: ArbFloat): ArbFloat;
  94. begin
  95. a0 := ((x + 1) * sin(x) + alfa * cos(x)) * spepow(x + 1, -alfa - 1);
  96. end;
  97. function ua0(u: ArbFloat): ArbFloat; {u=1/(x+1), of x=1/u-1}
  98. begin
  99. if u > 0 then
  100. ua0 := a0((1 - u) / u) / sqr(u)
  101. else
  102. ua0 := 0;
  103. end;
  104. function a1(x: ArbFloat): ArbFloat;
  105. var
  106. a: ArbFloat;
  107. begin
  108. a := spepow(x, -alfa);
  109. if scale then
  110. a1 := (alfa - 1) * a
  111. else
  112. a1 := a;
  113. end;
  114. function ua1(u: ArbFloat): ArbFloat; {u=ond/x of x=ond/u}
  115. begin
  116. if u > 0 then
  117. ua1 := a1(ond / u) * ond / sqr(u)
  118. else
  119. ua1 := 0;
  120. end;
  121. function a2(x: ArbFloat): ArbFloat;
  122. var
  123. a: ArbFloat;
  124. begin
  125. a := spepow(ln(x), -alfa) / x;
  126. if scale then
  127. a2 := (alfa - 1) * a
  128. else
  129. a2 := a;
  130. end;
  131. function ua2(u: ArbFloat): ArbFloat; {u=ond/x of x=ond/u}
  132. begin
  133. if u > 0 then
  134. ua2 := a2(ond / u) * ond / sqr(u)
  135. else
  136. ua2 := 0;
  137. end;
  138. function a3(x: ArbFloat): ArbFloat;
  139. var
  140. y: ArbFloat;
  141. begin
  142. if x = 0 then
  143. a3 := 0
  144. else
  145. begin
  146. y := spepow(x, alfa);
  147. a3 := alfa * y * sin(y) / (x * (x + 1)) + cos(y) / sqr(x + 1);
  148. end;
  149. end;
  150. function ua3(u: ArbFloat): ArbFloat; {u=1/(x+1), of x=1/u-1}
  151. begin
  152. if u > 0 then
  153. ua3 := a3((1 - u) / u) / sqr(u)
  154. else
  155. ua3 := 0;
  156. end;
  157. function a4(x: ArbFloat): ArbFloat;
  158. var
  159. y, z: ArbFloat;
  160. begin
  161. y := spepow(x, alfa);
  162. z := y * pi / 2;
  163. a4 := 2 * sin(z) / sqr(x + 1) - pi * alfa * y * cos(z) / (x * (x + 1));
  164. end;
  165. function ua4(u: ArbFloat): ArbFloat; {u=ond/x of x=ond/u}
  166. begin
  167. if u > 0 then
  168. ua4 := a4(ond / u) * ond / sqr(u)
  169. else
  170. ua4 := 0;
  171. end;
  172. function ss(x: ArbFloat): ArbFloat;
  173. var
  174. d, eps, r: ArbFloat;
  175. begin
  176. if x > 0.5 then
  177. begin
  178. d := frac(x);
  179. r := x - d;
  180. if d > 0.5 then
  181. begin
  182. d := 1 - d;
  183. r := r + 1;
  184. end;
  185. eps := 0.5 / spepow(r + 1, alfa);
  186. if d > eps then
  187. ss := 0
  188. else
  189. ss := (1 - d / eps) / (r * (r + 1) * eps);
  190. end
  191. else
  192. ss := 0;
  193. end;
  194. function uss(u: ArbFloat): ArbFloat; {u=ond/x of x=ond/u}
  195. begin
  196. if u > 0 then
  197. uss := ss(ond / u) * ond / sqr(u)
  198. else
  199. uss := 0;
  200. end;
  201. function sl(x: ArbFloat): ArbFloat;
  202. var
  203. sl1: ArbFloat;
  204. begin
  205. sl1 := sin(ln(x)) * spepow(x, -alfa);
  206. if scale then
  207. sl := sl1 / int1
  208. else
  209. sl := sl1;
  210. end;
  211. function usl(u: ArbFloat): ArbFloat; {u=ond/x of x=ond/u}
  212. begin
  213. if u > 0 then
  214. usl := sl(ond / u) * ond / sqr(u)
  215. else
  216. usl := 0;
  217. end;
  218. function se(x: ArbFloat): ArbFloat;
  219. var
  220. y, se1: ArbFloat;
  221. begin
  222. y := spepow(x, alfa);
  223. se1 := (sin(y) - alfa * (y / x) * cos(y)) * exp(-x);
  224. if scale then
  225. se := se1 / int1
  226. else
  227. se := se1;
  228. end;
  229. function use(u: ArbFloat): ArbFloat; {u=ond/x of x=ond/u}
  230. begin
  231. if u > 0 then
  232. use := se(ond / u) * ond / sqr(u)
  233. else
  234. use := 0;
  235. end;
  236. function v1(x: ArbFloat): ArbFloat;
  237. var
  238. a, y: ArbFloat;
  239. begin
  240. x := abs(x);
  241. alfa := abs(alfa);
  242. if x = 0 then
  243. begin
  244. if alfa = 1 then
  245. v1 := alfa / pi
  246. else
  247. v1 := 0;
  248. end
  249. else
  250. begin
  251. if x > 1 then
  252. a := -alfa - 1
  253. else
  254. a := alfa - 1;
  255. y := spepow(x, a);
  256. v1 := alfa * y / (pi * (sqr(x * y) + 1));
  257. end;
  258. end;
  259. function uv1(u: ArbFloat): ArbFloat; { u=«((2/ã)arctan(x)+1) of x=tan(«ã(2u-1)) }
  260. var
  261. y: ArbFloat; { 0 ó u ó 1 }
  262. begin
  263. if (u = 0) or (u = 1) then
  264. uv1 := 0
  265. else
  266. begin
  267. y := 1 / sqr(cos(pi * (u - 0.5)));
  268. uv1 := pi * v1(sqrt(y - 1)) * y;
  269. end;
  270. end;
  271. function v2(x: ArbFloat): ArbFloat;
  272. var
  273. v: ArbFloat;
  274. begin
  275. alfa := abs(alfa);
  276. if x > 0 then
  277. v := exp(-x / alfa)
  278. else
  279. if x < 0 then
  280. v := exp(x * alfa)
  281. else
  282. v := 1;
  283. if scale then
  284. v2 := v / (alfa + 1 / alfa)
  285. else
  286. v2 := v;
  287. end;
  288. function uv2(u: ArbFloat): ArbFloat; { u=«((2/ã)arctan(x)+1) of x=tan(«ã(2u-1)) }
  289. var
  290. y: ArbFloat; { 0 ó u ó 1 }
  291. begin
  292. if (u = 0) or (u = 1) then
  293. uv2 := 0
  294. else
  295. begin
  296. y := 1 / sqr(cos(pi * (u - 0.5)));
  297. if u > 0.5 then
  298. uv2 := pi * v2(sqrt(y - 1)) * y
  299. else
  300. uv2 := pi * v2(-sqrt(y - 1)) * y;
  301. end;
  302. end;
  303. var
  304. integral, ae, err: ArbFloat;
  305. term, num2: ArbInt;
  306. intex, First: boolean;
  307. procedure Header;
  308. var
  309. i: ArbInt;
  310. begin
  311. for i := 1 to 5 do
  312. if i = 3 then
  313. writeln(s: 3, ': ', Integraaltekst[u, i])
  314. else
  315. writeln('': 5, Integraaltekst[u, i]);
  316. end;
  317. procedure ShowResults;
  318. var
  319. f: ArbFloat;
  320. begin
  321. if First then
  322. writeln('alfa': num2, '': numdig - num2, 'ae': 7, ' ': 4, 'int': num2,
  323. '': numdig - num2, ' ', 'err': 7, ' ': 4, 'f': 6);
  324. First := False;
  325. if intex then
  326. f := inte - integral;
  327. case term of
  328. 1:
  329. begin
  330. Write(alfa: numdig, ae: 10, integral: numdig, ' ', err: 10, ' ');
  331. if intex then
  332. writeln(f: 10)
  333. else
  334. writeln;
  335. end;
  336. 2:
  337. begin
  338. Write(alfa: numdig, ae: 10, integral: numdig, ' ', err: 10, ' ');
  339. if intex then
  340. writeln(f: 10)
  341. else
  342. writeln;
  343. Writeln(' process afgebroken, te hoge nauwkeurigheid?');
  344. end;
  345. 3: Writeln('Verkeerde waarde ae (<=0) bij aanroep: ', ae: 8);
  346. 4:
  347. begin
  348. Write(alfa: numdig, ae: 10, integral: numdig, ' ', err: 10, ' ');
  349. if intex then
  350. writeln(f: 10)
  351. else
  352. writeln;
  353. writeln(' process afgebroken, moeilijk, mogelijk divergent?');
  354. end;
  355. end;
  356. end;
  357. const
  358. fint: array[boolean, 1..11] of rfunc1r =
  359. ((@ki, @a0, @a1, @a2, @a3, @a4, @ss, @sl, @se, @v1, @v2),
  360. (@uki, @ua0, @ua1, @ua2, @ua3, @ua4, @uss, @usl, @use, @uv1, @uv2));
  361. begin
  362. s := ParamStr(1);
  363. if s = '' then
  364. begin
  365. writeln(' Vergeten functienaam mee te geven!');
  366. writeln(' Kies uit: ', fnames);
  367. halt;
  368. end;
  369. for i := 1 to length(s) do
  370. s[i] := Upcase(s[i]);
  371. u := (Pos(s, fnames) + 2) div 3;
  372. if u = 0 then
  373. begin
  374. writeln(' Commandlineparameter ', s, ' bestaat niet');
  375. writeln(' Kies uit: ', fnames);
  376. halt;
  377. end;
  378. Write('program results int1fr function ' + s);
  379. case SizeOf(ArbFloat) of
  380. 4: writeln('(single)');
  381. 8: writeln('(double)');
  382. 6: writeln('(real)');
  383. end;
  384. num2 := numdig div 2;
  385. if Pos(s, 'a0 a4 a3 ss v1') > 0 then
  386. scale := True
  387. else
  388. begin
  389. Write(' scale ? (y or n)');
  390. readln(q);
  391. scale := Upcase(q) = 'Y';
  392. end;
  393. Write('Transformatie naar 0 => 1 ? (y or n)');
  394. readln(q);
  395. ond := ogs[u];
  396. f := fint[Upcase(q) = 'Y'][u];
  397. Header;
  398. Writeln('à en ae: ');
  399. First := True;
  400. while not eoln do
  401. begin
  402. Read(alfa, ae);
  403. intex := True;
  404. case u of
  405. 1: int1 := spebk0(alfa);
  406. 2:
  407. begin
  408. int1 := 1;
  409. intex := alfa > 0;
  410. end;
  411. 3:
  412. begin
  413. if alfa > 1 then
  414. int1 := 1 / (alfa - 1);
  415. intex := alfa > 1;
  416. end;
  417. 4:
  418. begin
  419. if alfa > 1 then
  420. int1 := 1 / (alfa - 1);
  421. intex := alfa > 1;
  422. end;
  423. 5:
  424. begin
  425. if alfa > 0 then
  426. int1 := 1
  427. else
  428. int1 := cos(1);
  429. intex := alfa > 0;
  430. end;
  431. 6:
  432. begin
  433. int1 := 1;
  434. intex := alfa > 0;
  435. end;
  436. 7: int1 := 1;
  437. 8:
  438. begin
  439. if alfa > 1 then
  440. int1 := 1 / (sqr(alfa - 1) + 1);
  441. intex := alfa > 1;
  442. end;
  443. 9: int1 := sin(1) / e;
  444. 10:
  445. begin
  446. int1 := 1;
  447. intex := alfa <> 0;
  448. end;
  449. 11:
  450. begin
  451. if alfa <> 0 then
  452. int1 := abs(alfa) + 1 / abs(alfa);
  453. intex := alfa <> 0;
  454. end;
  455. end;
  456. if scale then
  457. inte := 1
  458. else
  459. inte := int1;
  460. if Upcase(q) = 'Y' then
  461. int1fr(f, 0, 1, ae, integral, err, term)
  462. else if u < 10 then
  463. int1fr(f, ond, infinity, ae, integral, err, term)
  464. else
  465. int1fr(f, -infinity, infinity, ae, integral, err, term);
  466. Showresults;
  467. end;
  468. end.