SHOW:
|
|
- or go back to the newest paste.
1 | import java.io.File; | |
2 | import java.io.FileNotFoundException; | |
3 | import java.io.PrintWriter; | |
4 | ||
5 | import org.apache.commons.math.linear.Array2DRowRealMatrix; | |
6 | import org.apache.commons.math.linear.ArrayRealVector; | |
7 | import org.apache.commons.math.linear.DecompositionSolver; | |
8 | import org.apache.commons.math.linear.LUDecompositionImpl; | |
9 | import org.apache.commons.math.linear.RealMatrix; | |
10 | import org.apache.commons.math.linear.RealVector; | |
11 | ||
12 | ||
13 | public class Solution { | |
14 | public final static int a=-1; | |
15 | public final static int b=1; | |
16 | public final static int INTEGRAL_STEP=100000; | |
17 | public final static int lambda=1; | |
18 | private double h; | |
19 | private int n=1; | |
20 | private int sum; | |
21 | private double B[][]; | |
22 | private double d[]; | |
23 | private double c[]; | |
24 | private double eps_prev=0; | |
25 | public double f(double x){ | |
26 | return Math.cos(4*Math.PI*x)- | |
27 | (lambda*2*x*Math.sin(Math.PI*x))/(Math.PI*(x*x-16)); | |
28 | } | |
29 | public double K(double x, double s){ | |
30 | return Math.cos(Math.PI*x*s); | |
31 | //return Math.abs(x+s); | |
32 | } | |
33 | public double u(double x){ | |
34 | return Math.cos(4*Math.PI*x); | |
35 | //return Math.abs(x); | |
36 | } | |
37 | public double ksi_0(double x){ | |
38 | if (x>=a && x<=(a+h)){ | |
39 | return (a+h-x)/h; | |
40 | } | |
41 | else{ | |
42 | return 0; | |
43 | } | |
44 | } | |
45 | ||
46 | public double ksi_i(double x, int i){ | |
47 | if (x>=a+(i-1)*h && x<=(a+h*i)){ | |
48 | return (x-(a+(i-1)*h))/h; | |
49 | } | |
50 | else | |
51 | if (x>=a+i*h && x<=(a+h*(i+1))){ | |
52 | return (x+(i+1)*h-x)/h; | |
53 | } | |
54 | else{ | |
55 | return 0; | |
56 | } | |
57 | } | |
58 | ||
59 | public double ksi_n(double x){ | |
60 | if (x>=a+h*(n-1) && x<=(a+h*n)){ | |
61 | return (a+h-x)/h; | |
62 | } | |
63 | else{ | |
64 | return 0; | |
65 | } | |
66 | } | |
67 | public double intergrateBij(int i, int j){ | |
68 | double sum=0; | |
69 | double h=this.h/INTEGRAL_STEP;; | |
70 | double x; | |
71 | ||
72 | if (i==0){ | |
73 | for (int k=0; k<INTEGRAL_STEP;k++){ | |
74 | x=a+k*h; | |
75 | sum+=ksi_0(x)*K(x,a+this.h*j); | |
76 | } | |
77 | } | |
78 | else{ | |
79 | if (i==n){ | |
80 | for (int k=0; k<INTEGRAL_STEP;k++){ | |
81 | x=b-k*h; | |
82 | sum+=ksi_n(x)*K(x,a+this.h*j); | |
83 | } | |
84 | } | |
85 | else{ | |
86 | double x0=a+(i-1)*this.h; | |
87 | for (int k=0; k<INTEGRAL_STEP;k++){ | |
88 | x=x0+h*k; | |
89 | sum+=ksi_i(x,i)*K(x,a+this.h*j); | |
90 | } | |
91 | x0=a+i*this.h; | |
92 | for (int k=0; k<INTEGRAL_STEP;k++){ | |
93 | x=x0+h*k; | |
94 | sum+=ksi_i(x,i)*K(x,a+this.h*j); | |
95 | } | |
96 | } | |
97 | } | |
98 | return sum*h; | |
99 | } | |
100 | public double integrateDi(int i){ | |
101 | double sum=0; | |
102 | double h=this.h/INTEGRAL_STEP;; | |
103 | double x; | |
104 | if (i==0){ | |
105 | for (int k=0; k<INTEGRAL_STEP;k++){ | |
106 | x=a+k*h; | |
107 | sum+=ksi_0(x)*f(x); | |
108 | } | |
109 | } | |
110 | else{ | |
111 | if (i==n){ | |
112 | for (int k=0; k<INTEGRAL_STEP;k++){ | |
113 | x=b-k*h; | |
114 | sum+=ksi_n(x)*f(x); | |
115 | } | |
116 | } | |
117 | else{ | |
118 | double x0=a+(i-1)*this.h; | |
119 | for (int k=0; k<INTEGRAL_STEP;k++){ | |
120 | x=x0+h*k; | |
121 | sum+=ksi_i(x,i)*f(x); | |
122 | } | |
123 | x0=a+i*this.h; | |
124 | for (int k=0; k<INTEGRAL_STEP;k++){ | |
125 | x=x0+h*k; | |
126 | sum+=ksi_i(x,i)*f(x); | |
127 | } | |
128 | } | |
129 | } | |
130 | return sum*h; | |
131 | } | |
132 | public void initMartix(){ | |
133 | d=new double[n+1]; | |
134 | B=new double[n+1][]; | |
135 | for (int i=0; i<n+1; i++){ | |
136 | B[i]=new double[n+1]; | |
137 | } | |
138 | for (int i=0; i<n+1; i++){ | |
139 | for (int j=0; j<n+1; j++){ | |
140 | B[i][j]=lambda*intergrateBij(i, j); | |
141 | if (i!=j){ | |
142 | B[i][j]=-B[i][j]; | |
143 | } | |
144 | else{ | |
145 | B[i][j]=1-B[i][j]; | |
146 | } | |
147 | } | |
148 | d[i]=integrateDi(i); | |
149 | } | |
150 | RealMatrix coefficients=new Array2DRowRealMatrix(B); | |
151 | DecompositionSolver solver=new LUDecompositionImpl(coefficients).getSolver(); | |
152 | RealVector constants=new ArrayRealVector(d); | |
153 | RealVector solution=solver.solve(constants); | |
154 | c=solution.getData(); | |
155 | } | |
156 | public double u_wave(double x){ | |
157 | double sum=0; | |
158 | for (int i=0; i<n+1;i++){ | |
159 | sum+=c[i]*K(x,a+i*h); | |
160 | } | |
161 | sum*=lambda; | |
162 | sum+=f(x); | |
163 | return sum; | |
164 | } | |
165 | ||
166 | public void solve() throws FileNotFoundException{ | |
167 | PrintWriter writer; | |
168 | double x_i; | |
169 | double eps_max; | |
170 | double eps; | |
171 | double u_1,u_2; | |
172 | long start, finish; | |
173 | for (int i=1; i<=10; i++){ | |
174 | System.out.println("I="+i); | |
175 | start=System.currentTimeMillis(); | |
176 | sum=0; | |
177 | n*=2; | |
178 | h=(b-a)/(double)n; | |
179 | initMartix(); | |
180 | eps_max=-1; | |
181 | for (int j=1; j<=n; j++){ | |
182 | x_i=a+(j-1)*h; | |
183 | u_1=u(x_i); | |
184 | u_2=u_wave(x_i); | |
185 | System.out.println("with u="+u_1+" and u_wave="+u_2); | |
186 | eps=Math.abs(u_1-u_2); | |
187 | if (eps>eps_max){ | |
188 | eps_max=eps; | |
189 | } | |
190 | } | |
191 | finish=System.currentTimeMillis(); | |
192 | writer = new PrintWriter(new File("out" + i + ".txt")); | |
193 | writer.println("h_i = "+h); | |
194 | writer.println("eps = "+eps_max); | |
195 | writer.println("r_i = " + Math.log((eps_prev)/eps_max)/Math.log(2)); | |
196 | eps_prev=eps_max; | |
197 | writer.println("t_i = " + (finish-start)); | |
198 | writer.close(); | |
199 | } | |
200 | ||
201 | ||
202 | } | |
203 | ||
204 | } |