6-3 There is No Free Lunch (40分)

One day, CYLL found an interesting piece of commercial from newspaper: the Cyber-restaurant was offering a kind of “Lunch Special” which was said that one could “buy one get two for free”. That is, if you buy one of the dishes on their menu, denoted by d​i
​​ with price pi, you may get the two neighboring dishes di−1 and di+1 for free! If you pick up d1, then you may get d2and the last one d​n for free, and if you choose the last one d​n​​ , you may get d​n−1​​ and d​1​​ for free.

However, after investigation CYLL realized that there was no free lunch at all. The price p​i
​​ of the i-th dish was actually calculated by adding up twice the cost c​i
​​ of the dish and half of the costs of the two “free” dishes. Now given all the prices on the menu, you are asked to help CYLL find the cost of each of the dishes.

Format of function:
void Price( int n, double p[] );
where int n satisfies that 2< n ≤10000; double p[] is passed into the function as an array of prices p​,i and then will be returned storing the costs of the dishes.

Sample program of judge:

#include <stdio.h>

#define Max_size 10000 /* max number of dishes */

void Price( int n, double p[] );

int main()
{
    int n, i;
    double p[Max_size];

    scanf("%d", &n);
    for (i=0; i<n; i++) 
        scanf("%lf", &p[i]);
    Price(n, p);
    for (i=0; i<n; i++)
        printf("%.2f ", p[i]);
    printf("\n");

    return 0;
}

/* Your function will be put here */

Sample Input:
12 23.64 17.39 12.77 16.62 10.67 14.85 12.68 26.90 28.30 15.59 37.99 23.18
Sample Output:
9.20 5.58 3.24 7.00 1.99 6.36 2.25 10.01 11.52 0.50 17.65 4.88

尝试的方法：

高斯法：只能过前两个点，原因是内存不允许定义double * 10001 * 10001的数组

追赶法(?)：查了才知道有这个怪异的名字。主要想法就是$x_2$,…$x_{n-1}$都用$x_1$和$x_n$表示，然后带到最后一个方程里解。
只能过前两个点，原因是随着迭代次数增加，$x_1$,$x_n$的系数和常数项都会变的很大，超出double范围。（主要是因为这个矩阵的元素都是正的）

参数法+LU分解:
解法：
求解:$$\left[\begin{array}{ccccccc}2& 0.5& 0& ...& 0& 0& 0.5\\ & & & & & & \\ 0.5& 2& 0.5& ...& 0& 0& 0\\ & & & & & & \\ ...& ...& ...& ...& ...& ...& ...\\ & & & & & & \\ ...& ...& 0.5& 2& 0.5& ...& ...\\ & & & & & & \\ ...& ...& ...& ...& ...& ...& ...\\ & & & & & & \\ 0& 0& 0& ...& 0.5& 2& 0.5\\ & & & & & & \\ 0.5& 0& 0& ...& 0& 0.5& 2\end{array}\right]\left[\begin{array}{c}{x}_1\\ \\ x_2\\ \\ x_3\\ \\ ...\\ \\ ...\\ \\ ...\\ \\ x_n\end{array}\right]=\left[\begin{array}{c}{p}_1\\ \\ p_2\\ \\ p_3\\ \\ ...\\ \\ ...\\ \\ ...\\ \\ p_n\end{array}\right]$$
不妨设原矩阵为A,设左上方n-1阶矩阵为B,设$$c=[0.5,0,0,0...,0.5{]}^T$$
则有如下关系：
$$A=\left[\begin{array}{cc}B& c\\ c^T& 2\end{array}\right]$$
代入方程可得
$$\left[\begin{array}{c}B\left[\begin{array}{c}{x}_1\\ x_2\\ ...\\ x_{n-1}\end{array}\right]+c{x}_n\\ \\ c^T\left[\begin{array}{c}{x}_1\\ x_2\\ ...\\ x_{n-1}\end{array}\right]+2x_n\end{array}\right]=\left[\begin{array}{c}\left[\begin{array}{c}{p}_1\\ \\ p_2\\ \\ ...\\ \\ p_{n-1}\end{array}\right]\\ \\ p_n\end{array}\right]$$
取上方方块阵建立等式
不妨设$$x_i=u_i{x}_n+v_i,i=1,2,...,n$$代入等式
$$B(u{x}_n+v)+c{x}_n=P$$
$$Bv=P$$
$$Bu+c=0$$
就把原方程求解转化为解两个n-1阶三对角矩阵
$$\left[\begin{array}{ccccccc}2& 0.5& 0& ...& 0& 0& 0\\ & & & & & & \\ 0.5& 2& 0.5& ...& 0& 0& 0\\ & & & & & & \\ ...& ...& ...& ...& ...& ...& ...\\ & & & & & & \\ ...& ...& 0.5& 2& 0.5& ...& ...\\ & & & & & & \\ ...& ...& ...& ...& ...& ...& ...\\ & & & & & & \\ 0& 0& 0& ...& 0.5& 2& 0.5\\ & & & & & & \\ 0& 0& 0& ...& 0& 0.5& 2\end{array}\right]\left[\begin{array}{c}{u}_1\\ \\ u_2\\ \\ u_3\\ \\ ...\\ \\ ...\\ \\ ...\\ \\ u_{n-1}\end{array}\right]=\left[\begin{array}{c}-0.5\\ \\ 0\\ \\ ...\\ \\ 0\\ \\ ...\\ \\ 0\\ \\ -0.5\end{array}\right]$$
$$\left[\begin{array}{ccccccc}2& 0.5& 0& ...& 0& 0& 0\\ & & & & & & \\ 0.5& 2& 0.5& ...& 0& 0& 0\\ & & & & & & \\ ...& ...& ...& ...& ...& ...& ...\\ & & & & & & \\ ...& ...& 0.5& 2& 0.5& ...& ...\\ & & & & & & \\ ...& ...& ...& ...& ...& ...& ...\\ & & & & & & \\ 0& 0& 0& ...& 0.5& 2& 0.5\\ & & & & & & \\ 0& 0& 0& ...& 0& 0.5& 2\end{array}\right]\left[\begin{array}{c}{v}_1\\ \\ v_2\\ \\ v_3\\ \\ ...\\ \\ ...\\ \\ ...\\ \\ v_{n-1}\end{array}\right]=\left[\begin{array}{c}{p}_1\\ \\ p_2\\ \\ p_3\\ \\ ...\\ \\ ...\\ \\ ...\\ \\ p_{n-1}\end{array}\right]$$
最后利用
$$x_n=\frac{p_n-0.5v_1-0.5v_{n-1}}{0.5u_1+0.5u_{n-1}}$$解出$x_n$然后迭代求解。
解三对角矩阵：
用Court分解法

上代码

void CF(double* p, int n,double* x){
	double z[10001];
	double u[10001]; 
	z[0]=p[0]/2;
	u[0]=0.25;
	double l;
	for(int i=1;i<n-1;i++){
		l=2-0.5*u[i-1];
		u[i]=0.5/l;
		z[i]=(p[i]-0.5*z[i-1])/l;
	}
	l=2-0.5*u[n-2];
	z[n-1]=(p[n-1]-0.5*z[n-2])/l;
	x[n-1]=z[n-1];
	for(int i=n-2;i>=0;i--){
		x[i]=z[i]-u[i]*x[i+1];
	}
}
void Price( int n, double p[] ){
	double p1[10001];
	p1[0]=p1[n-2]=-0.5;
	for(int i=1;i<n-2;i++) p1[i]=0;
	double u[10001],v[10001];
	CF(p1,n-1,u);
	CF(p,n-1,v);
	double x[10001];
	x[n-1]=(p[n-1]-0.5*v[0]-0.5*v[n-2])/(2+0.5*u[0]+0.5*u[n-2]);
	for(int i=0;i<n-1;i++){
		x[i]=u[i]*x[n-1]+v[i];
	}
	for(int i=0;i<n;i++){
		p[i]=x[i];
	}
}

参考：

https://wenku.baidu.com/view/262c84ac941ea76e59fa042c.html#
https://blog.csdn.net/qq_32096491/article/details/78213914

---

祝大家解题愉快！