View difference between Paste ID: rfrmaS6K and iP62ACde
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
}