dbg := proc()
local HAM, TOL, x, N, h, M, t, alpha, f1, f2, f3, g1, g2, j, u1, u2, u3, v1, v2, i,
kf1, kf2, kf3, kg1, kg2, U;
1 HAM := [1];
2 Digits := 30;
3 TOL := 1/1000;
4 x[0] := 0;
5 N := 5000;
6 h := .1e-2;
7 M := 20;
8 t[0] := .6;
9 alpha := 0.;
10 f1 := (x, u1, u2, u3, v1, v2) -> u2;
11 f2 := (x, u1, u2, u3, v1, v2) -> u3;
12 f3 := (x, u1, u2, u3, v1, v2) -> u2^2-2*u1*u3-v1^2;
13 g1 := (x, u1, u2, v1, v2) -> v2;
14 g2 := (x, u1, u2, v1, v2) -> 2*u2*v1-2*u1*v2;
15 for j from 0 while j <= M do
16 u1[0] := 0;
17 u2[0] := t[j];
18 u3[0] := 0;
19 v1[0] := .8;
20 v2[0] := 0;
21 for i from 0 while i <= N-1 do
22 DEBUG(is((i = 221)*`and`(j = 2)));
23 kf1[i,1] := f1(x[i],u1[i],u2[i],u3[i],v1[i],v2[i]);
24 kf2[i,1] := f2(x[i],u1[i],u2[i],u3[i],v1[i],v2[i]);
25 kf3[i,1] := f3(x[i],u1[i],u2[i],u3[i],v1[i],v2[i]);
26 kg1[i,1] := g1(x[i],u1[i],u2[i],u3[i],v1[i],v2[i]);
27 kg2[i,1] := g2(x[i],u1[i],u2[i],u3[i],v1[i],v2[i]);
28 kf1[i,2] := f1(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,1],u2[i]+1/2*h*kf2[
i,1],u3[i]+1/2*h*kf3[i,1],v1[i]+1/2*h*kg1[i,1],v2[i]+1/2*h*
kg2[i,1]);
29 kf2[i,2] := f2(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,1],u2[i]+1/2*h*kf2[
i,1],u3[i]+1/2*h*kf3[i,1],v1[i]+1/2*h*kg1[i,1],v2[i]+1/2*h*
kg2[i,1]);
30 kf3[i,2] := f3(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,1],u2[i]+1/2*h*kf2[
i,1],u3[i]+1/2*h*kf3[i,1],v1[i]+1/2*h*kg1[i,1],v2[i]+1/2*h*
kg2[i,1]);
31 kg1[i,2] := g1(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,1],u2[i]+1/2*h*kf2[
i,1],u3[i]+1/2*h*kf3[i,1],v1[i]+1/2*h*kg1[i,1],v2[i]+1/2*h*
kg2[i,1]);
32 kg2[i,2] := g2(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,1],u2[i]+1/2*h*kf2[
i,1],u3[i]+1/2*h*kf3[i,1],v1[i]+1/2*h*kg1[i,1],v2[i]+1/2*h*
kg2[i,1]);
33 kf1[i,3] := f1(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,2],u2[i]+1/2*h*kf2[
i,2],u3[i]+1/2*h*kf3[i,2],v1[i]+1/2*h*kg1[i,2],v2[i]+1/2*h*
kg2[i,2]);
34 kf2[i,3] := f2(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,2],u2[i]+1/2*h*kf2[
i,2],u3[i]+1/2*h*kf3[i,2],v1[i]+1/2*h*kg1[i,2],v2[i]+1/2*h*
kg2[i,2]);
35 kf3[i,3] := f3(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,2],u2[i]+1/2*h*kf2[
i,2],u3[i]+1/2*h*kf3[i,2],v1[i]+1/2*h*kg1[i,2],v2[i]+1/2*h*
kg2[i,2]);
36 kg1[i,3] := g1(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,2],u2[i]+1/2*h*kf2[
i,2],u3[i]+1/2*h*kf3[i,2],v1[i]+1/2*h*kg1[i,2],v2[i]+1/2*h*
kg2[i,2]);
37 kg2[i,3] := g2(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,2],u2[i]+1/2*h*kf2[
i,2],u3[i]+1/2*h*kf3[i,2],v1[i]+1/2*h*kg1[i,2],v2[i]+1/2*h*
kg2[i,2]);
38 kf1[i,4] := f1(x[i]+h,h*kf1[i,3]+u1[i],h*kf2[i,3]+u2[i],h*kf3[i
,3]+u3[i],h*kg1[i,3]+v1[i],h*kg2[i,3]+v2[i]);
39 kf2[i,4] := f2(x[i]+h,h*kf1[i,3]+u1[i],h*kf2[i,3]+u2[i],h*kf3[i
,3]+u3[i],h*kg1[i,3]+v1[i],h*kg2[i,3]+v2[i]);
40 kf3[i,4] := f3(x[i]+h,h*kf1[i,3]+u1[i],h*kf2[i,3]+u2[i],h*kf3[i
,3]+u3[i],h*kg1[i,3]+v1[i],h*kg2[i,3]+v2[i]);
41 kg1[i,4] := g1(x[i]+h,h*kf1[i,3]+u1[i],h*kf2[i,3]+u2[i],h*kf3[i
,3]+u3[i],h*kg1[i,3]+v1[i],h*kg2[i,3]+v2[i]);
42 kg2[i,4] := g2(x[i]+h,h*kf1[i,3]+u1[i],h*kf2[i,3]+u2[i],h*kf3[i
,3]+u3[i],h*kg1[i,3]+v1[i],h*kg2[i,3]+v2[i]);
43 u1[i+1] := u1[i]+1/6*h*(kf1[i,1]+2*kf1[i,2]+2*kf1[i,3]+kf1[i,4]
);
44 u2[i+1] := u2[i]+1/6*h*(kf2[i,1]+2*kf2[i,2]+2*kf2[i,3]+kf2[i,4]
);
45 u3[i+1] := u3[i]+1/6*h*(kf3[i,1]+2*kf3[i,2]+2*kf3[i,3]+kf3[i,4]
);
46 v1[i+1] := v1[i]+1/6*h*(kg1[i,1]+2*kg1[i,2]+2*kg1[i,3]+kg1[i,4]
);
47 v2[i+1] := v2[i]+1/6*h*(kg2[i,1]+2*kg2[i,2]+2*kg2[i,3]+kg2[i,4]
);
48 if abs(v1[i+1]-1) <= TOL then
49 N := i+1;
50 print("Value of v = ",[(i+1)*h, v1[i+1], u2[N]]);
51 if abs(u2[i+1]) <= TOL then
52 print("Value of f' = ",u2[i+1]);
53 break
else
54 if j = 0 then
55 t[j+1] := t[j]+(alpha-u1[i+1])/(i+1)/h;
56 U[j] := u2[i+1];
57 print("Value of u1[j] = ",%)
else
58 U[j] := u2[i+1];
59 print("Value of u1[j] = ",%);
60 t[j+1] := t[j]-(U[j]-alpha)*(t[j]-t[j-1])/(U[j]-U[j
-1]);
61 print("Value of t = ",%)
end if
end if
end if
end do;
62 print(j,i)
end do
end proc
|
|