基于C#的三次插值

数值分析的上机作业，一开始想用Matlab做，但是发现方向不对，所以就换成了C#，其余的部分大都是窗体设计以及控件的编码，是个人而定，这里只展示下

核心代码

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
using System.Windows.Forms;

namespace Three_Times_spline_Interpolation
{
    class TTSIFuction
    {
    public float[] x={0}, y={0};
    public int edgeCondition;        //1为第一种边界条件，2为第一种边界条件，3为第一种边界条件
    public float head, tail;
    public int NodeNum;

    public void Cauc(Panel P)   //定义数组，图像编辑
    {
        float[] i = new float[NodeNum];
        float[] j = new float[NodeNum];
        float[] k = new float[NodeNum+1];        
        float[] i_j_k = new float[3];

        //for (int a = 0; a <= NodeNum; a++)
        //{
        //    x[a] = float.Parse(X[a]);
        //    y[a] = float.Parse(Y[a]);
        //}
        k[0] = ((y[1] - y[0]) / (x[1] - x[0]) - head) * 6 / (x[1] - x[0]); //头节点的6倍的差商，g（0）
        k[NodeNum] = (tail - (y[NodeNum] - y[NodeNum - 1]) / (x[NodeNum] - x[NodeNum - 1])) * 6 / (x[NodeNum] - x[NodeNum - 1]);    //g（Num）
        for (int b = 1; b< NodeNum; b++)
        {
            i[b] = (x[b] - x[b - 1]) / (x[b + 1] - x[b - 1]);   //求u
            j[b] = 1 - i[b];    //求人
            k[b] = ((y[b + 1] - y[b]) / (x[b + 1] - x[b]) - (y[b] - y[b - 1]) / (x[b] - x[b - 1])) * 6 / (x[b + 1] - x[b - 1]);
        }
        i_j_k[0] = (x[NodeNum] - x[NodeNum - 1]) / (x[1] - x[0] + x[NodeNum] - x[NodeNum - 1]);     //以下是三次样条求解
        i_j_k[1] = 1 - i_j_k[0];
        i_j_k[2] = (i_j_k[0] * (y[1] - y[0]) / (x[1] - x[0])+i_j_k[1]*(y[NodeNum]-y[NodeNum-1])/(x[NodeNum]-x[NodeNum-1]))*3;
        drawPic(x, y, SolveEquation(i, j, k, NodeNum, edgeCondition,i_j_k), P);

    }

    float[] SolveEquation(float[] i, float[] j, float[] k, int n, int edgeCondition, float[] i_j_k)
    {
        float[] M = new float[n + 1];

        for (int a = 0; a < n; a++)
        {
            M[a] = 0;
        }

        switch (edgeCondition)
        {
            case 1:         //第一种边界条件

                for(int p = 1;p<100;p++)
                {
                    //temp = M[0];
            
                    M[0] = (float)0.5 * (k[0] - M[1]);
                    for(int a=1;a<=n-1;a++)
                    {
                        M[a]= (float)0.5*(k[a]-i[a]*M[a-1]-j[a]*M[a+1]);
                    }
                    M[n]=(float)0.5*(k[n]-M[n-1]);
                    //e=M[0]-temp;
                    p++;
                }
                break;

            case 2:         //第二种边界条件

                k[1] = k[1] - i[0] * head;
                k[n - 1] = k[n - 1] - j[n - 1] * tail;
                i[1] = 0;
                j[n - 1] = 0;
                M[0] = head;
                M[n] = tail;
                for (int p = 1; p < 100; p++)
                {
                    
                    for (int a = 1; a <= n - 1; a++)
                    {
                        M[a] = (float)0.5 * (k[a] - i[a] * M[a - 1] - j[a] * M[a + 1]);
                    }
                    
                }

                break;
            case 3:         //第三种边界条件
                float[] Q = new float[n];
        
                for(int a =0;a<n;a++)
                {
                    Q[a]=0;
                }

                for (int p = 1; p < 100; p++)
                {
                    Q[0] = (float)0.5 * (k[1] - i[0] * Q[n - 1] - j[0] * Q[1]);
                    for (int a = 1; a <= n - 2; a++)
                    {
                        Q[a] = (float)0.5 * (k[a+1] - i[a] * Q[a - 1] - j[a] * Q[a + 1]);
                    }
                    Q[n - 1] = (float)0.5 * (i_j_k[2]-i_j_k[1]*Q[0]-i_j_k[0]*Q[n-2]);
                }
                break;
        }
        return M;
    }

    void drawPic(float[] X, float[] Y, float[] M, Panel P)      //画图了
    {
        try
        {
            float y;
            double y_temp=0;
            Graphics gr = P.CreateGraphics();
            Brush br1 = new SolidBrush(Color.Blue);
            Brush br2 = new SolidBrush(Color.Red);
            Pen pen1 = new Pen(br1, 11);
            Pen pen2 = new Pen(br2, 11);
            pen1.Width = (float)0.2;
            pen2.Width = (float)0.2;
            gr.DrawLine(pen1, 0,500, 1000,500);
            gr.DrawLine(pen1, 200, 0, 200, 1000);
            gr.DrawLine(pen1, 199, 250, 201, 250);
            gr.DrawLine(pen1, 250, 499, 250, 501);
            gr.DrawLine(pen1, 199, 450, 201, 450);
            gr.DrawLine(pen1, 700, 499, 700, 501);
            PointF pt_1 = new PointF(0, 0);
            PointF pt_2 = new PointF(0, 0);
            for (int q = 0; q <= NodeNum - 1; q++)
            {
                float temp_x = X[q], temp_y = Y[q];
                for (float x = (float)(X[q] + 0.0005); x <= X[q + 1]; x = (float)(x + 0.0005))
                {
             //       y = (float)((M[q] * Math.Pow((X[q + 1] - x), 3) / 6 + M[q + 1] * Math.Pow((x - X[q]), 3) / 6 + (Y[q] - M[q] * Math.Pow((X[q + 1] - X[q]), 2) / 6) * (X[q + 1] - x) + (Y[q + 1] - M[q + 1] * Math.Pow((X[q + 1] - X[q]), 2) / 6) * (x - X[q])) / (X[q + 1] - X[q]));
                    y_temp = M[q] * Math.Pow((X[q + 1] - x), 3) ;
                    y_temp = y_temp + M[q + 1] * Math.Pow((x - X[q]), 3) ;
                    y_temp = y_temp + (6*Y[q] - M[q] * Math.Pow((X[q + 1] - X[q]), 2) ) * (X[q + 1] - x);
                    y_temp = y_temp + (6*Y[q + 1] - M[q + 1] * Math.Pow((X[q + 1] - X[q]), 2) ) * (x - X[q]);
                    y = (float)((y_temp / (X[q + 1] - X[q]))/6);
                    pt_1.X = temp_x*50+200 ;
                    pt_1.Y = temp_y*(-50)+500 ;
                    pt_2.X = x*50+200 ;
                    pt_2.Y = y*(-50)+500 ;

                    gr.DrawLine(pen2, pt_1, pt_2);
                    temp_x = x;
                    temp_y = y;
                }
            }
        }
        catch { Exception ex; }
    }
    }
}