1
2

#归一化 收到0~1之间
#正则化不是归一化 

1
2
3
4
5
6
7
8

from sklearn.preprocessing import MinMaxScaler
 
data = [[-1, 2], [-0.5, 6], [0, 10], [1, 18]]
 
#不太熟悉numpy的小伙伴，能够判断data的结构吗？
#如果换成表是什么样子？
import pandas as pd
pd.DataFrame(data)

1
2
3

from sklearn import tree
from sklearn.datasets import load_wine
from sklearn.model_selection import train_test_split

1

wine = load_wine()

1

wine.data.shape

1

wine.target

1
2

import pandas as pd
pd.concat([pd.DataFrame(wine.data),pd.DataFrame(wine.target)],axis=1)

1
2
3
4
5

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split     #分割训练集和测试集
from sklearn.neighbors import KNeighborsClassifier       #K近邻

1
2
3

iris=datasets.load_iris()
iris_X = iris.data
iris_y=iris.target

1

iris_X[:2,:]#四个属性两朵花

1

iris_y#三类花

1
2

#分割训练集和测试集
X_train,X_test,y_train,y_test=train_test_split(iris_X,iris_y,test_size=0.3)#测试占30%

1

y_train#顺便打乱了数据

1
2

knn = KNeighborsClassifier() #使用K近邻分类
knn.fit(X_train,y_train)     #输入测试集

1

knn.predict(X_test)#预测

1

y_test

1
2

from sklearn import datasets #导入基础数据
from sklearn.linear_model import LinearRegression  #导入模型

1
2
3
4

#导入数据
loaded_data = datasets.load_boston()
data_X = loaded_data.data
data_y = loaded_data.target

1
2
3

#定义模型并训练
model = LinearRegression()#线性模型
model.fit(data_X,data_y)

1
2

#预测
model.predict(data_X[:4,:])

1

data_y[:4]

1
2
3

import numpy as np
from sklearn.tree import DecisionTreeRegressor #回归树
import matplotlib.pyplot as plt

1
2
3

#生成噪声曲线
rng = np.random.RandomState(1)
rng.rand(2,3)#生成区分行列的二维数据，一维数据不分行列，生成0~1之间的数

1
2
3

X = np.sort(5*rng.rand(80,1),axis=0)
y = np.sin(X).ravel()#导入二维数据，ravel降维y，可将任意维度降维为一维
plt.scatter(X,y,s=20,edgecolors="black",c="darkorange",label="data")

1
2

y[::5] += 3 * (0.5 - rng.rand(16))#每5个数字加随机+-0.5之间的数当成噪声
plt.scatter(X,y,s=20,edgecolors="black",c="darkorange",label="data")

1
2
3
4
5

#设定参数的模型导入并训练
regr_1 = DecisionTreeRegressor(max_depth=2)
regr_2 = DecisionTreeRegressor(max_depth=5)
regr_1.fit(X,y)
regr_2.fit(X,y)

1
2

X_test = np.arange(0.0,5.0,0.01)[:,np.newaxis]
#取0~5以0.01为步长的有顺序点，后面将1维升维2，跟reshape（-1，1）作用一致

1
2

y_1 = regr_1.predict(X_test)
y_2 = regr_2.predict(X_test)

1
2
3
4
5
6
7
8
9

plt.figure()
plt.scatter(X,y,s=20,edgecolors="black",c="darkorange",label="data")#边框，点，图例名称
plt.plot(X_test,y_1,color="cornflowerblue",label="max_depth=2",linewidth=2)
plt.plot(X_test,y_2,color="yellowgreen",label="max_depth=5",linewidth=2)
plt.xlabel("data")
plt.ylabel("target")
plt.title("Decision Tree Regression")
plt.legend()#显示图例
plt.show()

1
2

import os
os.sys.executable

1

'/usr/bin/python3'

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

import os
libs = {
    "requests","jieba","beautifulsoup4","matplotlib","numpy","pandas","openpyxl","tensorflow","仿照格式这里填写要安装的包"
}
try:
    for lib in libs:
        #这里的地址是上一步显示的解释器位置
        os.system('/usr/bin/python3 -m pip install '+lib)
        print("Successful")
except:
    print('error')

1
2
3
4
5

import numpy as np
import pandas as pd


data = pd.read_excel(r'C:\Users\祈LHL\Desktop\data.xlsx')

1

data.head()

1

y = data.sort_values(by="y-coordinate" , ascending=False)

1
2
3
4
5
6

step = 4.892/6
i = 0
y_list = []
section = []
#y6  x3

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

#分y
while i < 6:
    step_y = float("%.8f"%(i*step))
    a_y = step_y-2.446
    y_min = y.loc[y['y-coordinate'] >= a_y]
    i +=1
    step_y = float("%.8f"%(i*step))
    b_y = step_y-2.446
    y_max = y_min.loc[y_min['y-coordinate'] < b_y]
    y_list.append(y_max)
    section_y = (a_y,b_y)
    section.append(section_y)
#section[5]

1
2
3
4
5

j = 0
step = 1.562/3
x_list = []
section_x_list = []
section_y_list = []

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

#在y中分x
while j <= 5:
    temp = y_list[j]
    m = 0
    while m <= 2:
        step_x = float("%.8f"%(m*step))
        a_x = step_x-12.639
        x_min = temp.loc[temp['x-coordinate'] >= a_x]
        m +=1
        step_x = float("%.8f"%(m*step))
        b_x = step_x-12.639
        x_max = x_min.loc[temp['x-coordinate'] < b_x]
        x_list.append(x_max)
        section_x = (a_x,b_x)
        section_x_list.append(section_x)
        section_y_list.append(section[j])
    j += 1
#18个块

1
2
3

section_list = []
Q_list = []
i = 0

1
2
3
4
5

while i <= 17:
    temp_l = x_list[i]
    Q = (temp_l["density"]*temp_l["z-velocity"]*temp_l["z-face-area"]).sum()
    i +=1
    Q_list.append(Q)

1

section_x_list

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

[(-12.639, -12.118333329999999),
 (-12.118333329999999, -11.597666669999999),
 (-11.597666669999999, -11.077),
 (-12.639, -12.118333329999999),
 (-12.118333329999999, -11.597666669999999),
 (-11.597666669999999, -11.077),
 (-12.639, -12.118333329999999),
 (-12.118333329999999, -11.597666669999999),
 (-11.597666669999999, -11.077),
 (-12.639, -12.118333329999999),
 (-12.118333329999999, -11.597666669999999),
 (-11.597666669999999, -11.077),
 (-12.639, -12.118333329999999),
 (-12.118333329999999, -11.597666669999999),
 (-11.597666669999999, -11.077),
 (-12.639, -12.118333329999999),
 (-12.118333329999999, -11.597666669999999),
 (-11.597666669999999, -11.077)]

1

section_y_list

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

[(-2.446, -1.63066667),
 (-2.446, -1.63066667),
 (-2.446, -1.63066667),
 (-1.63066667, -0.8153333300000001),
 (-1.63066667, -0.8153333300000001),
 (-1.63066667, -0.8153333300000001),
 (-0.8153333300000001, 0.0),
 (-0.8153333300000001, 0.0),
 (-0.8153333300000001, 0.0),
 (0.0, 0.8153333299999996),
 (0.0, 0.8153333299999996),
 (0.0, 0.8153333299999996),
 (0.8153333299999996, 1.6306666699999997),
 (0.8153333299999996, 1.6306666699999997),
 (0.8153333299999996, 1.6306666699999997),
 (1.6306666699999997, 2.446),
 (1.6306666699999997, 2.446),
 (1.6306666699999997, 2.446)]

1
2
3
4
5
6

index = []
i = 0
while i <= 2:
    index.append(section_x_list[i])
    i+=1
index

1
2
3

[(-12.639, -12.118333329999999),
 (-12.118333329999999, -11.597666669999999),
 (-11.597666669999999, -11.077)]

1
2
3
4
5
6

cols = []
i = 0
while i <= 17:
    cols.append(section_y_list[i])
    i+=3
cols

1
2
3
4
5
6

[(-2.446, -1.63066667),
 (-1.63066667, -0.8153333300000001),
 (-0.8153333300000001, 0.0),
 (0.0, 0.8153333299999996),
 (0.8153333299999996, 1.6306666699999997),
 (1.6306666699999997, 2.446)]

1
2
3
4
5
6
7

j = 0
data = {}
while j <=5:
    col = cols[j]
    data[col]= 1
    j += 1
data

1
2
3
4
5
6

{(-2.446, -1.63066667): 1,
 (-1.63066667, -0.8153333300000001): 1,
 (-0.8153333300000001, 0.0): 1,
 (0.0, 0.8153333299999996): 1,
 (0.8153333299999996, 1.6306666699999997): 1,
 (1.6306666699999997, 2.446): 1}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

temp_list = Q_list.copy()
r = 0
l = 0
for i in data:
    content = []
    for j in range(3):
        content.append(Q_list[l])
        l+=1
    data[i] = content
data

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

{(-2.446, -1.63066667): [12.427282344087539,
  12.86438295520224,
  8.104732674789329],
 (-1.63066667, -0.8153333300000001): [14.059159859218617,
  14.0170078290342,
  6.50613023081204],
 (-0.8153333300000001, 0.0): [13.851341644545903,
  14.102576555137208,
  6.549793869809068],
 (0.0, 0.8153333299999996): [13.873905539737082,
  14.100577509379482,
  6.539991338980014],
 (0.8153333299999996, 1.6306666699999997): [14.030126072372612,
  13.966730831575635,
  6.467075694188603],
 (1.6306666699999997, 2.446): [11.644887356127308,
  12.004108151505823,
  7.422577844039829]}

1
2

df = pd.DataFrame(data=data, index=index)
df

  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
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114

%% SIR差分模型
frame=importdata('美国covid19疫情数据l-history改.csv');
date_row=size(frame.data,1);
data=frame.data;
date=frame.textdata(:,1);
real_date=date(2:date_row+1,:);
E=zeros(2,90);
l=0.000125;    %日接触率
m=0.01;   %日治愈率0.0005
E(1,1)=data(345,12)/300000000;
E(2,1)=1-E(1,1);
for i=1:90
    E(1,i+1)=l*E(2,i)-m*E(1,i)+E(1,i);
    E(2,i+1)=E(2,i)-l*E(1,i)*E(2,i);
end
X=flip(data(311:345,12)');% 5.31  4.2 
rate=flip(data(311:345,12)')%现阳性率
rate(isnan(rate))=0;
n=length(rate);
rt=0:1:n-1;
%a0=[100,10]; %初值
figure(1)
h1 = plot(rt,X/300000000,'*'); %画点
hold on;
a=E(1,:);
b=E(2,:);
h2 = plot(a,'r')
xlabel('日期');         %设置横坐标名
ylabel('感染人数');   %设置纵坐标名
legend([h1 h2],'3~5月感染人数','SIR模型迭代曲线','Location','NorthWest');
hold on
%%  预测
X=flip(data(255:345,12)');%311行后为第一题
rate=flip(data(255:345,12)')%现阳性率
rate(isnan(rate))=0;
n=length(rate);
rt=0:1:n-1;
figure(2)
h3 = plot(rt,X/300000000,'*'); %画点
hold on;
a=E(1,:);
b=E(2,:);
h4=plot(a,'r')
xlabel('日期');
ylabel('感染人数'); 
legend([h3 h4],'3~6月感染人数','SIR模型迭代曲线','Location','NorthWest');
hold on
%% 接种疫苗前
frame=importdata('usc.csv');
date_row=size(frame.data,1);
data=frame.data;
date=frame.textdata(:,1);
real_date=date(2:date_row+1,:);
E=zeros(3,150);
l=0.0018;    %日接触率 0.0018
m=0.01;   %日治愈率 0.01
E(1,1)=data(320,1)/300000000;
E(3,1)=28664448/300000000%43714928  33714928  28664448
E(2,1)=1-E(1,1)-E(3,1);

for i=1:150
    E(1,i+1)=E(1,i)+l*E(2,i)-m*E(1,i);
    E(2,i+1)=E(2,i)-l*E(1,i)*E(2,i);
    E(3,i+1)=1-E(1,i+1)- E(2,i+1);
end
X=data(320:486,1)';%80-482  3.1-7.1
rate=data(320:486,1)'%
rate(isnan(rate))=0;
n=length(rate);
rt=0:1:n-1;
figure(6)
h11 = plot(rt,X/300000000,'*'); %画点
hold on;
a=E(1,:);
b=E(2,:);
c=E(3,:);
t=1:151;
h1=plot(t,a,'r')%感染人数
%plot(t,b,'b')%健康人数
plot(t,c,'r')%移除者
%legend('i(t)','s(t)','r(t)')
hold on
%% 接种疫苗后
E=zeros(3,150);
l=0.0018;    %日接触率 0.0018
m=0.01;   %日治愈率 0.01
E(1,1)=data(320,1)/300000000;
E(3,1)=53714928/300000000%43714928  33714928  28664448
E(2,1)=1-E(1,1)-E(3,1);
for i=1:150
    E(1,i+1)=E(1,i)+l*E(2,i)-m*E(1,i);
    E(2,i+1)=E(2,i)-l*E(1,i)*E(2,i);
    E(3,i+1)=1-E(1,i+1)- E(2,i+1);
end
X=data(320:486,1)';%80-482  3.1-7.1
rate=data(320:486,1)'%
rate(isnan(rate))=0;
n=length(rate);
rt=0:1:n-1;
figure(6)
%h11 = plot(rt,X/300000000,'*'); %画点
hold on;
a=E(1,:);
b=E(2,:);
c=E(3,:);
t=1:151;
h2=plot(t,a,'b')%感染人数
%plot(t,b,'b')%健康人数
plot(t,c,'b')%移除者
%legend('i(t)','s(t)','r(t)')
hold on
xlabel('日期');         %设置横坐标名
ylabel('r(t)与i（t）');   %设置纵坐标名
legend([h1 h2 h11],'接种疫苗前','接种疫苗后','感染人数原始数据','Location','NorthWest');

 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

%改进的指数增长模型
x=flip([324985536 322941312 320635168 318300992 315993728 313830976 311556864 309321664 306771520 304093952 301231200 298379904 295516608 292805312 290107936 287625184 284968960 282162400 279040000 275854016 272656992])
y=flip([0.006330017,0.007192424,0.007333235,0.007301613,0.006891455,0.007299188,0.007226135,0.008312845,0.008805068,0.009503504,0.009555925,0.00968912,0.009259723,0.009297836,0.008631901,0.009321099,0.009946612,0.011189794,0.011549529,0.011725443,0.012112401])
n=length(x);
t=0:1:n-1;
rk=zeros(1,n);
rk(1)=(-3*x(1)+4*x(2)-x(3))/2;
rk(n)=(x(n-2)-4*x(n-1)+3*x(n))/2;
for i=2:n-1
rk(i)=(x(i+1)-x(i-1))/2;
end
rk=rk./x;
p=polyfit(t,rk,2);
r0=p(1);
r1=p(2);
r2=p(3);
x0=x(1);
R=r0*t.^2+r1*t+r2;
X=x0*exp((r0*t.^3)/3+(r1*t.^2)/2+r2*t);
figure(1)
hold on;
xlabel('year');         %设置横坐标名
ylabel('rate');         %设置纵坐标名
grid on                 %网格线
plot(t,y,'r*')
figure(2)
hold on;
xlabel('year');         %设置横坐标名
ylabel('population');   %设置纵坐标名
grid on                 %网格线
plot(t,x,'r*',t,X)
figure(3)
hold on;
xlabel('year');         %设置横坐标名
ylabel('rate');         %设置纵坐标名
grid on                 %网格线
plot(t,y,'r*',t,R')
%2018
t1=t(21)+1;
x2018=x0*exp((r0*t1.^3)/3+(r1*t1.^2)/2+r2*t1)
y2018=r0*t1.^2+r1*t1+r2
%2019
t1=t(21)+2;
x2019=x0*exp((r0*t1.^3)/3+(r1*t1.^2)/2+r2*t1)
y2019=r0*t1.^2+r1*t1+r2

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

clc;
clear all;
syms x t;
f(t)=6*(exp(-0.1155*t)-exp(-0.1386*t));
h(t)=exp(-0.1386*t);
figure(1);
fplot(f,[0 25]);
hold on;
fplot(h,[0,25]);
grid on;
xlabel('时间t/h');
ylabel('药量');
legend('血液中的药量y/mg','胃肠道中的药量x/mg');
figure(2)
fplot(f,[0,25]);
hold on;
fplot(h,[0,25]);
f(t)=-(exp(-0.693*t)-exp(-0.1386*t))/4;
fplot(f,[0,25]);
xlabel('时间t/h');
ylabel('药量');
legend('血液中的药量y/mg','胃肠道中的药量x/mg','施救后血液中的药量y/mg');

  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
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

clear;
close all;
clc;
pho=[300;862;74.2;1.18];
c=[1377;2100;1726;1005];
lamda=[0.082;0.37;0.045;0.028];
a=lamda./(pho.*c);
 
d=[0.6;6;3.6;5]*10^-3;
 
TT=273.15;
T_in=37;
T_out=75;
T_s=48.08;
 
xmin=0;
xmax=sum(d);
 
N=5400;
h=0.05*10^-3;
k=1;
r=k/h^2;
I=round((xmax-xmin)/h);
 
A=zeros(1,I);
B=zeros(1,I+1);
C=zeros(1,I);
 
N1=round(d(1)/h);
N2=round(d(2)/h);
N3=round(d(3)/h);
N4=round(d(4)/h);
 
for i=1:N1
    A(i)=-a(1)*r;
    B(i)=2+2*r*a(1);
    C(i)=-r*a(1);
end
for i=N1+1:N1+N2
    A(i)=-a(2)*r;
    B(i)=2+2*r*a(2);
    C(i)=-r*a(2);
end
for i=N1+N2+1:N1+N2+N3
    A(i)=-a(3)*r;
    B(i)=2+2*r*a(3);
    C(i)=-r*a(3);
end
for i=N1+N2+N3+1:N1+N2+N3+N4
    A(i)=-a(4)*r;
    B(i)=2+2*r*a(4);
    C(i)=-r*a(4);
end
 
T=zeros(I+1,N+1);
T(:,1)=(T_in+TT)*ones(I+1,1);
 
T_xt=xlsread('CUMCM-2018-Problem-A-Chinese-Appendix.xlsx');
 
h_min=110;
h_max=120;
delta_h=0.1;
H1=h_min:delta_h:h_max;
delta=zeros(1,length(H1));
 
for j=1:length(H1)
    h1=h_min+(j-1)*delta_h;
    k1=lamda(1);k2=lamda(2);k3=lamda(3);k4=lamda(4);
    x1=d(1);x2=d(1)+d(2);x3=d(1)+d(2)+d(3);x4=d(1)+d(2)+d(2)+d(4);
    t1=T_out+TT;t2=T_in+TT;t3=T_s+TT;
    
    h5=-((h1*k2*k3*k4*t1)/(k1*k2*k3*k4-h1*k1*k2*k3*x3-h1*k1*k2*k4*x2 ...
        -h1*k1*k3*k4*x1+h1*k1*k2*k3*x4+h1*k1*k2*k4*x3+h1*k1*k3*k4*x2+h1*k2*k3*k4*x1)-(h1*k2*k3*k4*t3)...
        /(k1*k2*k3*k4-h1*k1*k2*k3*x3-h1*k1*k2*k4*x2-h1*k1*k3*k4*x1+h1*k1*k2*k3*x4+h1*k1*k2*k4*x3+h1*k1*k3*k4*x2+h1*k2*k3*k4*x1))...
        /(t2/k1-t3/k1);
    
    AA=diag(B)+diag(A,1)+diag(C,-1);
    AA(1,1)=lamda(1)/h+h1;
    AA(1,2)=-lamda(1)/h;
    AA(I+1,I)=-lamda(4)/h;
    AA(I+1,I+1)=lamda(4)/h+h5;
    
    AA(N1+1,N1)=-lamda(1);
    AA(N1+1,N1+1)=lamda(1)+lamda(2);
    AA(N1+1,N1+2)=-lamda(2);
    
    AA(N1+N2+1,N1+N2)=-lamda(2);
    AA(N1+N2+1,N1+N2+1)=lamda(2)+lamda(3);
    AA(N1+N2+1,N1+N2+2)=-lamda(3);
    
    AA(N1+N2+N3+1,N1+N2+N3)=-lamda(3);
    AA(N1+N2+N3+1,N1+N2+N3+1)=lamda(3)+lamda(4);
    AA(N1+N2+N3+1,N1+N2+N3+2)=-lamda(4);
    
    for n=1:k:N
        D=zeros(I+1,1);
        D(1)=h1*(T_out+TT);
        D(I+1)=h5*(T_in+TT);
        for i=2:1:N1
            D(i)=r*a(1)*T(i-1,n)+(2-2*r*a(1))*T(i,n)+r*a(1)*T(i+1,n);
        end
        for i=N1+1:1:N1+N2
            D(i)=r*a(2)*T(i-1,n)+(2-2*r*a(2))*T(i,n)+r*a(2)*T(i+1,n);
        end
        for i=N1+N2+1:1:N1+N2+N3
            D(i)=r*a(3)*T(i-1,n)+(2-2*r*a(3))*T(i,n)+r*a(3)*T(i+1,n);
        end
        for i=N1+N2+N3+1:1:N1+N2+N3+N4
            D(i)=r*a(4)*T(i-1,n)+(2-2*r*a(4))*T(i,n)+r*a(4)*T(i+1,n);
        end
        D(N1+1)=0;
        D(N1+N2+1)=0;
        D(N1+N2+N3+1)=0;
        T(:,n+1)=AA\D;
    end
   delta(j)=sqrt(sum((T_xt(:,2)-T(end,:)'+TT).^2)/length(T_xt(:,1)));
end
%图二 
figure(1);
mesh(0:k:N,1000*(0:h:sum(d)),(T-TT));
%图三
T_problem1=zeros(N+1,4);
T_problem1(:,1)=T(1,:)';
T_problem1(:,2)=T(N1+1,:)';
T_problem1(:,3)=T(N2+N1+1,:)';
T_problem1(:,4)=T(N3+N2+N1+1,:)';
T_problem1=T_problem1-TT;
figure(2);
plot(0:k:N,T_problem1(:,1)',0:k:N,T_problem1(:,2)',0:k:N,T_problem1(:,3)',0:k:N,T_problem1(:,4)',0:k:N,T_xt(:,2)');