最小二乘
This question will develop a set of functions to least square fit the linear model 𝑦=𝑘𝑥+𝑞 to arbitrary data provided in an input file, i.e. identify the coefficients 𝑘 and 𝑞 to optimally overlap the data points (𝑥, 𝑦) available in the input file.
本问题将开发一套函数，对输入文件中的任意数据进行线性模型𝑦=𝑘𝑥+𝑞的最小二乘拟合，即识别出系数𝑘和𝑞，使输入文件中可获得的数据点(𝑥，𝑦)最优重叠。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
typedef struct {
char *filename;
int num_points;
double *x;
double *y;
} data;
typedef struct {
double k;
double q;
double r2;
} fitcoef;
/* 4-4 */
int data_size(char filename[]){
FILE *fp = NULL;
int ans=0;
double x,y;
fp = fopen(filename, "r");
if(fp == NULL){
return -1;
}
while(fscanf(fp,"%lf %lf",&x,&y)!=EOF){
ans++;
}
fclose(fp);
return ans;
}
/* 4-5 */
int data_read(data* m_data){
FILE *fp = NULL;
int idx=0;
m_data->x=(double*) malloc((m_data->num_points)*sizeof(double));
m_data->y=(double*) malloc((m_data->num_points)*sizeof(double));
fp = fopen(m_data->filename, "r");
if(fp == NULL){
return -1;
}
for(int i=0;i<m_data->num_points;i++){
fscanf(fp,"%lf %lf",&(m_data->x[idx]),&(m_data->y[idx]));
idx++;
}
fclose(fp);
return 0;
}
/* 4-6 */
fitcoef linear_fit(data* m_data){
fitcoef ans;
double _x = 0.0, _y = 0.0;
double s_xx = 0.0, s_yy = 0.0, s_yx = 0.0;
for(int i=0;i<m_data->num_points;i++){
_x+=m_data->x[i];
_y+=m_data->y[i];
}
_x/=m_data->num_points;
_y/=m_data->num_points;
for(int i=0;i<m_data->num_points;i++){
s_xx+=pow(m_data->x[i] - _x,2.0);
}
s_xx/=m_data->num_points;
for (int i = 0; i < m_data->num_points; i++) {
s_yy += pow((m_data->y[i] - _y),2.0) ;
}
s_yy/=m_data->num_points;
for (int i = 0; i < m_data->num_points; i++) {
s_yx += (m_data->x[i] - _x) * (m_data->y[i] - _y);
}
s_yx/=m_data->num_points;
ans.k=s_yx/s_xx;
ans.q=_y-ans.k*_x;
ans.r2=s_yx*s_yx/s_xx/s_yy;
return ans;
}
/* 4-7 */
void print_data(data* m_data){
printf("\nData file: %s\n", m_data->filename);
printf("Number of data points: %m_data\n\n", m_data->num_points );
printf("Data points (x_i, y_i):\n");
for(int i=0;i<m_data->num_points;i++){
printf("i=%d:\t (%lf, %lf)\n", i,m_data->x[i],m_data->y[i]);
}
}
int main()
{
data points;
fitcoef fc;
points.filename = "numbers.txt";
points.num_points = data_size(points.filename);
data_read(&points);
print_data(&points);
printf("\nFitting line y = kx + q to the data...\n");
fc = linear_fit(&points);
printf("\nFit coeficients:\n");
printf("k = %lf\n", fc.k);
printf("q = %lf\n", fc.q);
printf("\nCorrelation coefficient:\n");
printf("R2 = %lf\n\n",fc.r2);
return 0;
}
- Create an input file named numbers.txt containing the following two columns of numbers:
创建一个名为numbers.txt的输入文件，包含以下两列数字:
0 0.8
1 2.2
2 2.9
3 3.8
4 3.9
where the numbers in each row are separated by white space.
每一行的数字用空格隔开。 - Implement the structure data as follows:
如下所示, 实现结构data: - The members k and q will hold the linear coefficients obtained from the fitting, and r2 is the so-called linear correlation coefficient describing the quality of the fit. Good fit leads r2 close to 1, r2 is always smaller than or equal to 1.
成员k和q将保留拟合得到的线性系数，r2是所谓的描述拟合质量的线性相关系数。好的拟合使r2接近1,r2总是小于等于1。 - Write a function data_size(…) which takes as an input a two-column data filename such as our numbers.txt and returns an integer corresponding to the total number of data points (𝑥,𝑦) in the data file. Assume that a general data file with arbitrary number of points (rows) can be passed to the function and usenumbers.txt as a test case when your function data_size should return number 5.
编写一个函数data_size(…)，它接受一个两列数据文件名(如numbers.txt)作为输入，并返回一个与数据文件中数据点总数(𝑥，𝑦)相对应的整数。 假设一个具有任意数量点(行)的通用数据文件可以作为一个测试用例传递给函数和usenumbers.txt，而函数data_size应该返回5。 - Write a function data_read(…) which takes as input a pointer to structure data. Assume that the members filename and num_points of the input structure have already been assignment values (see the function main example below). The function data_read should then populate the structure member arrays x and y with the data available in the file filename. The new structure - now also storing data points (𝑥, 𝑦) - should be returned through the input pointer to data. The function data_read also returns a type int equal either to 0 if the file open/close operations were successful or −1 if they were unsuccessful.
编写一个函数data_read(…)，将指向结构体data的指针作为输入。假设输入结构体filename和num_points的成员已经赋值(参见下面函数main的例子)。 函数data_read应该用文件filename中可用的数据填充结构成员数组x和y。新的结构——现在也存储数据点(𝑥，𝑦)——应该通过指向data的输入指针返回。 函数data_read还返回类型int，如果文件打开/关闭操作成功，则返回0;如果文件打开/关闭操作失败，则返回−1。 - To perform linear fit of 𝑦=𝑘𝑥+𝑞 to the data (𝑥, 𝑦), write a function linear_fit(…) which takes as input a pointer to structure data and returns structure of type fitcoef. The function should internally evaluate the following expressions:
为了对数据(𝑥，𝑦)进行𝑦=𝑘𝑥+𝑞的线性拟合，编写一个函数linear_fit(…)，以结构data的指针作为输入，返回类型fitcoef的结构。该函数应在内部对以下表达式求值:
- Write a function print_data(…) which returns type void. It takes as an input a pointer to structure data. The function should print the content of the structure data as shown in the example below in part 9). Hint: the skeleton of function print_data is as follows:
编写一个返回类型void的函数print_data(…)，它接受一个指向结构data的指针作为输入。该函数应该打印结构data的内容，如下面第9部分中的示例所示。提示:函数print_data的框架如下:
向量旋转
In this question we will study applications of rotation matrices to rotate vectors by a specified angle.
在这个问题中，我们将研究旋转矩阵的应用，使向量旋转一个指定的角度。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define M_PI 3.14159265358979323846
#define M_PI_2 (M_PI/2.0)
#define N 2
/* 5-1 */
double vector_norm(double* v){
return sqrt(pow(v[0],2.0)+ pow(v[1], 2.0));
}
/* 5-2 */
void unit_vector(double* v){
double l=vector_norm(v);
v[0]/=l;v[1]/=l;
}
/* 5-3 */
void rotation_matrix(double** R, double theta){
R[0][0]=cos(theta);R[1][0]=sin(theta);R[0][1]=-R[1][0];R[1][1]=R[0][0];
}
/* 5-4 */
void rotate_vector(double** R,double v[2]){
double tmp_v[2];
tmp_v[0] = v[0];
tmp_v[1] = v[1];
v[0]= tmp_v[0]*R[0][0]+ tmp_v[1]*R[1][0];
v[1]= tmp_v[0]*R[0][1]+ tmp_v[1]*R[1][1];
}
/* 5-5 */
double rotation_sequence(double* v,double theta,int n){
double d_t=-theta/n;
double** r =(double **)malloc(N*sizeof(double*));
for (int i=0;i<N;i++){
r[i]=(double *)malloc(N*sizeof(double));
}
rotation_matrix(r,d_t);
double ans=0.0;
double last[2],now[2],d[2];
last[0]=v[0];
last[1]=v[1];
for(int i=0;i<n;i++){
now[0]=last[0];
now[1]=last[1];
rotate_vector(r,now);
d[0]=now[0]-last[0];
d[1]=now[1]-last[1];
last[0]=now[0];
last[1]=now[1];
ans += vector_norm(d);
}
v[0]=last[0];
v[1]=last[1];
return ans;
}
/* 5-6 */
double** find_rotation_matrix(double* v,double* w){
double t=(v[0]*w[0]+v[1]*w[1])/(pow(vector_norm(v),2.0));
t=acos(t);
double** ans =(double **)malloc(2*sizeof(double*));
for (int i=0;i<2;i++){
ans[i]=(double *)malloc(2*sizeof(double));
}
rotation_matrix(ans,t);
return ans;
}
/* 5-7 */
void print_vector(double *v) {
printf("[%5.3lf, %5.3lf]\n", v[0], v[1]);
}
void print_matrix(double **R) {
printf("|%5.3lf\t%5.3lf|\n", R[0][0], R[0][1]);
printf("|%5.3lf\t%5.3lf|\n", R[1][0], R[1][1]);
}
int main()
{
double** R =(double **)malloc(N*sizeof(double*));
for (int i=0;i<N;i++){
R[i]=(double *)malloc(N*sizeof(double));
}
double u[N] = {1.0, 0.0};
double v[N] = {1.0, 1.0};
double w[N] = {0.0, 1.0};
double theta = M_PI_2/2.0; // M_PI_2 from math.h
int num_theta = 100;
double arc_len;
printf("\nVector u:\n");
print_vector(u);
printf("\nRotation angle: %lf rad (%5.2lf deg)", theta, theta*180.0/M_PI);
rotation_matrix(R, theta);
printf("\nRotation matrix R:\n");
print_matrix(R);
rotate_vector(R, u);
printf("\nVector u rotated by R:\n");
print_vector(u);
printf("\nVector v:\n");
print_vector(v);
printf("Vector w:\n");
print_vector(w);
printf("\nMatrix rotating v to w:\n");
R = find_rotation_matrix(v, w);
print_matrix(R);
printf("\n");
v[0]=0.707;
v[1]=0.707;
printf("Rotating vector v: ");
print_vector(v);
arc_len = rotation_sequence(v, theta, num_theta);
printf("by %lf rad (%5.2lf deg)", theta,theta*180.0/M_PI);
printf(" in %d steps.\n", num_theta);
printf("Traversed arc length: %lf.\n", arc_len);
printf("The final vector is: ");
print_vector(v);
printf("\n");
return 0;
}
- Write a function vector_norm(…) which returns type double and takes as input an array v of two elements of type double. The function should compute and return the norm (length) |v| of the vector v.
编写一个函数vector_norm(…)，返回类型double，并接受包含两个类型double的元素的数组v作为输入。该函数应该计算并返回向量v的范数(长度) |v|。 - Write a function unit_vector(…) which returns type void and takes as an input an array v of two elements of type double. The function should normalize the vector v to unit length, and return the normalized unit vector through the input array v. You can use the function vector_norm developed in part 1) to solve this question.
编写一个函数unit_vector(…)，返回类型void，并接受包含两个double类型元素的数组v作为输入。 函数应该将向量v归一化为单位长度，并通过输入数组v返回归一化的单位向量。您可以使用在第1部分中开发的函数vector_norm来解决这个问题。 - Write a function rotation_matrix(…) which returns type void and takes two inputs. The first input is a two-dimensional array R having 2𝑥2 elements of type double. The second input is a variable theta of type double. The function rotation_matrix should compute the rotation matrix:
编写一个函数rotation_matrix(…)，返回类型void并接受两个输入。 第一个输入是一个二维数组R，有2个类型为double的2𝑥2元素。 第二个输入是一个类型为double的变量theta。 函数rotation_matrix应该计算旋转矩阵:
where 𝜃 is set equal to the input variable theta and represents the rotation angle in radians. The function should return the rotation matrix corresponding to theta through the input argument R.
其中𝜃被设为输入变量theta，表示以弧度表示的旋转角度。 函数应该通过输入参数R.返回与theta对应的旋转矩阵。
- Write a function rotate_vector(…) which returns type void and takes two inputs. The first input is a two-dimensional array R having 2𝑥2 elements of type double representing the rotation matrix as determined in part 3) above. The second input is a two-element array v of type double representing a vector. The function rotate_vector should calculate a new rotated vector by computing the product of the matrix R and input vector v. The rotated vector should be returned through the input argument v.
编写一个函数rotate_vector(…)，返回void类型并接受两个输入。第一个输入是一个二维数组R，其中有2个类型为double的2𝑥2元素，表示上面第3部分中确定的旋转矩阵。 第二个输入是一个双元素数组v，类型为double，表示一个向量。 函数rotate_vector应该通过计算矩阵R和输入向量v的乘积来计算一个新的旋转后的向量。旋转后的向量应该通过输入参数v返回。 - Write a function rotation_sequence(…) which returns a type double and takes three input arguments. The first input argument is a two-element array v of type double representing the vector to be rotated. The second argument is the maximum angle theta in radians by which the vector v will be rotated. The third argument is an integer n determining the number of incremental rotations towards the angle theta. For example, if theta = 𝜋/2 and 𝑛=4, we will require 4 sub-rotations by an angle 𝜋/8 to reach the full angle theta = 𝜋/2. The function rotation_sequence should return the distance L travelled by the tip of vector v during the n-step rotation towards theta. Note that in the limit of large n this distance should approach the value L = theta ∙|v|. You can use the functions implemented in previous parts to facilitate the solution. Hint: the distance L can be calculated by summing the vector differences v’-v, where v and v’ are the vectors before and after an incremental rotation.
编写一个函数rotation_sequence(…)，它返回一个double类型并接受三个输入参数。 第一个输入参数是一个双元素数组v，类型为double，表示要旋转的向量。 第二个参数是以弧度表示的旋转矢量v的最大角度theta。 第三个参数是一个整数n，它决定了朝向角度的增量旋转次数。 例如，如果theta =𝜋/2和𝑛=4，我们将需要4个子旋转一个角度𝜋/8，以达到theta =𝜋/2的全角度。 函数rotation_sequence应该返回矢量v的尖端在向着的n步旋转过程中所走过的距离L。 注意，在n较大的极限下，这一距离应该接近L = theta ∙|v|的值。 您可以使用前面部分中实现的函数来促进解决方案的实现。 提示:距离L可以通过矢量差v’-v求和得到，其中v和v’是增量旋转前后的矢量。 - Write a function find_rotation_matrix(…) which takes as input two arrays v and w each having two elements of type double. These arrays represent two vectors. The function should return a pointer-to-pointer to type double. This pointer is to point to an array identified within the function find_rotation_matrix representing the rotation matrix required to rotate vector v to w. You can assume that vectors v and w are such that v can be rotated towards w in anticlockwise direction, i.e. the rotation angle is positive. Again, previously developed functions can aid in solving this question. Hint: the rotation angle between unit vectors □(→┬a ) and □(→┬b ) can be found by using the dot product as cos𝜃=□(→┬a∙→┬b ).
编写一个函数find_rotation_matrix(…)，它接受两个数组v和w作为输入，每个数组都有两个double类型的元素。 这些数组表示两个向量。 该函数应返回一个指向double类型的指针到指针。 这个指针指向函数find_rotation_matrix中标识的数组，该数组表示将向量v旋转到w所需的旋转矩阵。 你可以假设向量v和w可以使v逆时针向w方向旋转, 即旋转角是正的。 同样，以前开发的函数可以帮助解决这个问题。 提示:单位向量□(→┬a )与□(→┬b )之间的旋转角度可用点积cos𝜃=□(→┬a∙→┬b )求得。
欧拉方法求洛伦兹方程
This question will solve the famous Lorenz equations to generate the strange attractor (see Fig. 1(b)) as a signature of chaotic system. Lorenz equations read:
这个问题将求解著名的洛伦兹方程，生成奇异吸引子(见图1(b))作为混沌系统的特征。洛伦茨方程:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define Ti 0.0 // initial time
#define Tf 50.0 // final time
#define dT 0.001 // time step
#define Y01 0.0 // initial condition 1
#define Y02 1.0 // initial condition 2
#define Y03 0.0 // initial condition 3
#define sigma 10.0 // Lorenz equations parameters
#define b 8.0/3.0
#define r 28.0
/* 6-1 */
void lorenz(double* args, double* vars) {
vars[0] = sigma * (args[1] - args[0]);
vars[1] = r * args[0] - args[1] - args[0] * args[2];
vars[2] = args[0] * args[1] - b * args[2];
}
/* 6-2 */
void euler(double** p, void (*f)(double*, double*), int n) {
double vars[3], args[3];
for (int i = 0; i < n; i++) {
args[0] = p[0][i];
args[1] = p[1][i];
args[2] = p[2][i];
f(args, vars);
p[0][i + 1] = p[0][i] + (vars[0] * dT);
p[1][i + 1] = p[1][i] + (vars[1] * dT);
p[2][i + 1] = p[2][i] + (vars[2] * dT);
}
}
/* 6-3 */
int main()
{
double* t;
double** y;
int N = (Tf - Ti) / dT;
t = (double*)malloc(N * sizeof(double));
y = (double**)malloc(N * sizeof(double*));
for (int i = 0; i < 3; i++) {
y[i] = (double*)malloc(N * sizeof(double));
}
y[0][0] = Y01;
y[1][0] = Y02;
y[2][0] = Y03;
euler(y, lorenz, N); // y to store solutions
for (int i = 0; i <= N; i++) {
t[i] = Ti + dT * i; // define time
}
FILE* fp =fopen("data.txt", "w");
for (int i = 0; i <N; i++) {
fprintf(fp, "%f\t%f\t%f\t%f\n", t[i], y[0][i], y[1][i], y[2][i]);
}
close(fp);
return 0;
}
- Write a function lorenz(…) which returns type void and takes as input two arrays of type double. The first array holds the values of variables 𝑥, 𝑦, 𝑧, and the second array holds the values of functions 𝑓, 𝑔, ℎ evaluated at those values. Thus the first array serves as input and second array serves as output of the function lorenz. The parameters 𝜎, 𝑟, and 𝑏 can be defined as symbolic constants (see the main program below).
编写一个函数lorenz(…)，返回类型void并接受两个类型double的数组作为输入。 第一个数组包含变量𝑥、𝑦、𝑧的值，第二个数组保存根据这些值计算的函数𝑓、𝑔、ℎ的值。 因此，第一个数组作为lorenz(…)函数的输入，第二个数组作为洛伦兹函数的输出。参数𝜎、𝑟和𝑏可以定义为符号常量(参见下面的main程序)。 - Write a function euler(…) which implements the Euler method to solve Lorenz equations defined above. The prototype of the function euler is given as:
写一个函数euler(…)实现了欧拉方法，以解决洛伦兹方程的上述定义。 函数euler的原型为:
void euler(double **p, double (*f)(double, double), int n)
The double-pointer ∗∗𝑝 holds an array of size 3×𝑛 to store the 𝑥, 𝑦, 𝑧 solutions over 𝑛 time instants, starting from initial time 𝑇𝑖 to final time 𝑇𝑓 and assuming time step 𝑑𝑇. The value of integer 𝑛 is passed to the function euler as the last argument. The second argument is expected to be the pointer to function lorenz implemented in part 1) above. Note that the function euler does not require passing the time array as input. This is because Lorenz equations do not contain explicit time dependence. The initial conditions and the time parameters can be defined as symbolic constants (see the main program below).
双指针∗∗𝑝保存一个大小为3×𝑛的数组，用于存储𝑛时间瞬间内的𝑥、𝑦、𝑧解，从初始时间𝑇𝑖到最终时间𝑇𝑓，假设时间步长𝑑𝑇。 整型𝑛的值作为最后一个参数传递给函数euler。第二个参数预计是上面第1部分中实现的lorenz函数的指针。注意，函数euler不需要传递时间数组作为输入。这是因为洛伦兹方程不包含显式的时间依赖关系。初始条件和时间参数可以定义为符号常量(参见下面的main程序)。
2265
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
