using MathNet Numerics LinearAlgebra Double using System using System

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
using MathNet.Numerics.LinearAlgebra.Double;
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Diagnostics;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
namespace hord_approximation {
public partial class Form1 : Form {
public Form1() {
InitializeComponent();
}
private void Form1_Paint(object sender, PaintEventArgs e) {
var g = e.Graphics;
var D = new double[,]{{50,200},{100,100},{150,100},{200,200}};
int k = 3;
var X = new int[] { 0, 0, 0, 1, 2, 2, 2 };
var N = new double[D.GetLength(0), 4];
var T = new double[] { 0, 0.817256, 1.182744, 2 };
for (var i = 0; i < D.GetLength(0); i++) {
for (var j = 0; j < 4; j++) {
N[i, j] = GetBasisFunc(j + 1, k, T[i], X);
}
}
var Dmatrix = DenseMatrix.OfArray(D);
var Nmatrix = DenseMatrix.OfArray(N);
var NInverse = Nmatrix.Inverse();
var Bmatrix = NInverse * Dmatrix;
var B = Bmatrix.ToArray();
for (var i = 0; i < B.GetLength(0); i++) {
g.DrawEllipse(new Pen(Color.Blue), (float)B[i, 0], (float)B[i, 1], 4f, 4f);
}
for (double t = 0; t <= 2; t += 0.001) {
double x = B[0, 0] * GetBasisFunc(1, k, t, X) + B[1, 0] * GetBasisFunc(2, k, t, X) + B[2, 0] * GetBasisFunc(3, k, t, X) + B[3, 0] * GetBasisFunc(4, k, t, X);
double y = B[0, 1] * GetBasisFunc(1, k, t, X) + B[1, 1] * GetBasisFunc(2, k, t, X) + B[2, 1] * GetBasisFunc(3, k, t, X) + B[3, 1] * GetBasisFunc(4, k, t, X);
g.DrawEllipse(new Pen(Color.Red), (float)x, (float)y, 2f, 2f);
}
}
private double GetBasisFunc(int i, int k, double t, int[] X) {
double result = 0;
if (k == 1) {
if ((t >= X[(i) - 1]) && (t <= X[(i + 1) - 1]))
result = 1;
else
result = 0;
}
else {
if (X[(i + k - 1) - 1] - X[(i) - 1] != 0)
result = (t - X[(i) - 1]) * GetBasisFunc(i, k - 1, t, X) / (X[(i + k - 1) - 1] - X[(i) - 1]);
if (X[(i + k) - 1] - X[(i + 1) - 1] != 0)
result += (X[(i + k) - 1] - t) * GetBasisFunc(i + 1, k - 1, t, X) / (X[(i + k) - 1] - X[(i + 1) - 1]);
}
return result;
}
}
}