Differential-Equations

本文最后更新于：2023年6月19日 晚上

非刚性常微分方程

形如类的方程,反求出 y 的值,并最终画出原函数图像

ode23 – 求解非刚性微分方程 - 低阶方法

语法

[t,y] = ode23(odefun,tspan,y0)
[t,y] = ode23(odefun,tspan,y0,options)
[t,y,te,ye,ie] = ode23(odefun,tspan,y0,options)
sol = ode23(___)

说明

[t,y] = ode23(odefun,tspan,y0)（其中 tspan = [t0 tf]）求微分方程组 从 到 的积分，初始条件为 。解数组 中的每一行都与列向量 中返回的值相对应。,必须返回矩阵才能一一对应!
参数odefun
要求解的函数，指定为指向待积分函数的句柄。
对于标量 t 和列向量 y 来说，函数 dydt = odefun(t,y) 必须返回数据类型为 single 或 double 的列向量 dydt，该列向量对应于 f(t,y)。odefun 必须同时接受输入参数 t 和 y，即使其中一个参数未在函数中使用也是如此。
例如，要解算 y′=5y−3，请使用此函数：

1
2
function dydt = odefun(t,y)
dydt = 5*y-3;

对于方程组，odefun 的输出为向量。向量中的每个元素是一个方程的解。例如，要求解

使用函数：

1
2
3
4
function dydt = odefun(t,y)
dydt = zeros(2,1);
dydt(1) = y(1)+2*y(2);
dydt(2) = 3*y(1)+2*y(2);

tspan

积分区间，指定为向量。tspan 必须至少是一个二元素向量 [t0 tf]，用于指定初始时间和最终时间。要获取 t0 到 tf 之间的特定时间的解，请使用 [t0,t1,t2,…,tf] 形式的长向量。tspan 中的元素必须单调递增或单调递减。
求解器在初始时间 tspan(1) 施加由 y0 给出的初始条件，然后求 tspan(1) 到 tspan(end) 的积分：
如果 tspan 有两个元素，[t0 tf]，求解器将返回在该区间内的每个内部积分步计算的解。
如果 tspan 包含两个以上的元素，[t0,t1,t2,…,tf]，求解器将返回在给定点处计算的解。但是，求解器不会精确步进到 tspan 中指定的每个点。此时，求解器使用自己的内部积分步来计算解，然后在 tspan 中请求的各点处计算解。在指定点处生成的解与在每个内部积分步计算的解具有相同的准确度级别。示例

对于二阶非刚性 ODE 方程

如$y_1’’- (1-y_1^2)*y_1’+y_1=0$,这是一个van der Pol 方程.
二阶方程做一阶变换
令$y_2=y_1’$,因此$y_2’=(1-y_1^2)*y_2+y_1$
类似的对应关系
具体 matlab 实现

1
2
3
%vdp1.m文件
function dydt = vdp1(t,y)
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];%将方程两边作为二维矩阵的两个行向量

使用 ode23 函数、时间区间 [0 20] 和初始值 [2 ; 0]来解算该 ODE。生成的输出即为时间点 t 的列向量和解数组 y。y 中的每一行都与 t 的相应行中返回的时间相对应。y 的第一列与 相对应，第二列与 相对应。

虽然初始值是上下排列,但最后得到的数据是左右排列(一列是对应一个 y,感觉好像转置了)

调用文件

1
2
3
4
5
6
[t,y] = ode23(@vdp1,[0 20],[2; 0]);%使用@文件名可以调用里面的函数(最好文件名和函数名一致)
plot(t,y(:,1),'-o',t,y(:,2),'-o')%绘制y1和y2的解图
title('Solution of van der Pol Equation (\mu = 1) with ODE23');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

注意
ode23 仅适用于使用两个输入参数（t 和 y）的函数。但是，通过在函数外部定义参数并在指定函数句柄时传递这些参数，可以传入额外参数。

传入额外参数

解$y’’=\frac{A}{B}ty$
重写令$y_1’=y_2$,因此$y_2’=\frac{A}{B}ty_1$
函数文件

1
2
function dydt = odefun(t,y,a,b)
dydt = [y(2);(a/b)*t.*y(1)]

调用文件

1
2
3
4
5
6
A = 1;
B = 2;
tspan = [0 5];
y0 = [0;0.01];
[t,y] = ode23(@(t,y) odefun(t,y,A,B), tspan, y0);
plot(t,y(:,1),'-o',t,y(:,2),'-.')

或者一个文件试了一下也是可以的

1
2
3
4
5
6
7
odefun = @(t,y,a,b)[y(2);(a/b)*t.*y(1)];
A = 1;
B = 2;
tspan = [0 5];
y0 = [0;0.01];
[t,y] = ode23(@(t,y) odefun(t,y,A,B), tspan, y0);
plot(t,y(:,1),'-o',t,y(:,2),'-.')

总结

ode23 是 Bogacki 和 Shampine 的显式 Runge-Kutta (2,3) 对的实现。在容差较宽松且刚度适中的情况下，它可能比 ode45 更加有效。ode23 是单步求解器

ode45 – 求解非刚性微分方程 - 中阶方法

语法

[t,y] = ode45(odefun,tspan,y0)
[t,y] = ode45(odefun,tspan,y0,options)
[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)
sol = ode45(___)
ode45 是一个通用型 ODE 求解器，是解算大多数问题时的首选。但是，对于刚性问题或需要较高准确性的问题，其他 ODE 求解器可能更适合。
其实用法和 ode23 几乎一样,只是生成的数据图的数据点更密了.

刚性常微分方程

区别

所谓刚性方程，就是说存在两(多)重尺度，一个尺度比另外一个尺度大很多。所导致的麻烦就是在计算中很难兼顾两者。例如下面的方程：dx=-100x-100.1ydy=100.1x-100y 两个特征值 lambda_1 = -200.1，lambda_2 = -0.1 所以解表现为 a1_exp(-200.1_X) + a2_exp(-0.1_X)，无论你用什么样的尺度(单一尺度)都不能很好刻画解的行为。一个是快变行为，一个是慢变行为。 所有这样的方程计算时候，稳定性条件比较苛刻。实际情况要比这还复杂得多。
对于刚性和非刚性微分方程的区分，可以简单的转变为在将原方程转换为常微分方程组后，进行一个简单的系数判断：
例如：

$$
y’’’ - 3y’’ - y’y = 0, y(0) = 0, y’(0) = 0, y’’(0) = -1;
$$

在这里可以设  y1 = y, y2 = y’, y3 = y’’, 有

$$
 y_1’ = y_2,      　y_1(0) = 0,
$$

$$
 y_2’ = y_3,       y_2(0) = 1,
$$

$$
 y_3’ = 3y_3 + y_2y_1,   y_3(0) = -1,
$$　　
这里可以简单判断出方程组的右侧系数矩阵值差异不大，得到的特征值差异随之也不大，可以简单判断为非刚性微分方程。
MATLAB中解非刚性微分方程常用 ode45 ;

又例如：

在这里可以设, 有


　　
这里可以简单判断出方程组的右侧系数矩阵值差异较大，得到的特征值差异随之较大，可以简单判断为刚性微分方程。
MATLAB中解刚性微分方程常用 ode15s , ode23s , ode23t , ode23tb ;

如果区分不出来就多试试.

1
2
dy=@(t,y)[y(2);y(3);3*y(3)+y(2)*y(1)];
[T,Y]=ode23s(@(t,y)dy(t,y),[0 1],[0;1;-1]);%函数传参一定要声明句柄,参数前后都要写

解常微分方程符号解

常微分方程在 Matlab 中按如下规定重新表达： 符号 D 表示对变量的求导。Dy 表示对变量 y 求一阶导数，当需要求变量的 n 阶导 数时，用 Dn 表示，D4y 表示对变量 y 求 4 阶导数。由此，常微分方程 $y’’+2y’= y$ 在 Matlab 中，将写成

求解常微分方程的通解

无初边值条件的常微分方程的解就是该方程的通解。
dsolve(' diff_equation'，'var')
式中 diff_equation 为待解的常微分方程，第 1 种格式将以变量 t 为自变量进行求解， 第 2 种格式则需定义自变量 var。
栗子:
解
程序如下:

1
2
3
4
5
6
syms x y
diff_equ='x^2+y+(x-2*y)*Dy=0';%字符串形式的变量!注意!
dsolve(diff_equ,'x')
%ans =
% x/2 + ((4*x^3)/3 + x^2 + C4)^(1/2)/2
% x/2 - ((4*x^3)/3 + x^2 + C4)^(1/2)/2

求解常微分方程的初边值问题

dsolve('diff_equation'，'condition1，condition2，…'，'var')其中 condition1，condition2，… 即为微分方程的初边值条件。
例如y=dsolve('D3y-D2y=x','y(1)=8,Dy(1)=7,D2y(2)=4','x')
结果为

1
2
y =
x*((exp(-1)*(19*exp(1) - 14))/2 - 1) + 7*exp(-2)*exp(x) - x^2/2 - x^3/6 + (exp(-1)*(19*exp(1) - 14))/2 - (exp(-1)*(25*exp(1) - 21))/3 - 1

求解常微分方程组

``dsolve``(``'diff_equ1，diff_equ2，…'``，``'``var``'``)``%用于求通解
``dsolve``(``'diff_equ1，diff_equ2，…'``，``'``condition1，condition2，…``'``，``'``var``'``)``%用于求初始值
例如,求

这个方程组的通解和在初始值为的解.

1
2
3
4
equ1='D2f+3*g=sin(x)';
equ2='Dg+Df=cos(x)';
[general_f,general_g]=dsolve(equ1,equ2,'x')
[f,g]=dsolve(equ1,equ2,'Df(2)=0,f(3)=3,g(5)=1','x')

目前没有覆盖到多变量的常微分方程数值解,也不知道ode45传初始值的矩阵究竟用列向量还是行向量(似乎都可以,无论逗号或是分号)

实战学习

双变量初始值下的常微分方程组

伏尔泰拉(Volterra)弱肉强食模型

问题提出

问题解决

即初值问题

其中,

  • r是食饵独立生存的时候自然增长率

  • a 是由于捕食者存在导致减少的比例系数

  • c是捕食者独自存在时候的死亡率

  • b是由于食饵存在导致死亡率减少的比例系数

它们均大于0.

注意多变量的矩阵变换
或许必须将x’,y’放在方程组的一边用于迭代

但是如果x’与y’是乘的关系呢?
似乎不用考虑,因为高数里面也没有涉及到
函数文件

1
2
3
4
5
6
7
8
function dd = shier(t,x)
r = 1;
c =0.5;
a = 0.1;
b = 0.02;
dd = diag([r-a*x(2,:),-c+b*x(1,:)])*x;%把原来的x,y放在一个矩阵里就可以只用一个符号变量进行索引
%对于左边是x'的,就把x初始值乘进去,然后索引y的初始值,其实如果反过来索引也可以
dd返回的是[x',y']',然后交给ode45处理

调用文件

1
2
3
tspan = 0:0.1:15;
[t,x] = ode45(@shier,tspan,[25;2]);%我认为t只是传进来迭代次数,与常微分方程组本身并无关系
plot(t,x);

结果

这里甚至可以看出图像具有周期性,那么怎么计算出周期呢?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
tspan = 0:0.1:15;
[t,x] = ode45(@shier,tspan,[25;2]);%我认为t只是传进来迭代次数,与常微分方程组本身并无关系
%方法一
%a = find((abs(x(1,1)-x(1:length(x),1)))<1);%如果是min的话他只会索引到第一个原值
这里
%方法二
a = find(x(:,1)==(x(1,1)-min(abs((x(1,1)-x(2:length(x),1))))));%这里不是加就是减,一定会索引到最小值的,多试两遍就出来了

figure(1)
%plot(t,x,t(a(2,1),1),x(a(2,1),1),'*');%方法一对应的找数据点方法
plot(t,x,t(a,1),x(a,1),'*');%方法二对应的找数据点的方法
hold on;
line([0 15],[x(a,1) x(a,1)],'LineStyle','--');%画一条水平线
line([t(a,1) t(a,1)],[0 x(a,1)],'LineStyle','--');%画一条竖直线
figure(2);
plot(x(:,1),x(:,2))

效果图


其中dd参数在传入初始值[25;2]后运算结果(即t=0.1时)

1
2
3
4
5
ans =

 25*r - 25*a*y
 2*b*x - 2*c
 %这是一个列向量,有两行,其中y为2,x为25,这里只是为了便于理解
  而结果最终返回的x是一个n*2的矩阵,第一列是x,第二列是y

结果的数值表示为

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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
ans =

                   0  25.000000000000000   2.000000000000000
   0.100000000000000  27.081808053754724   2.004112695948303
   0.200000000000000  29.334409623386446   2.016970388233990
   0.300000000000000  31.768915238930827   2.039429632085728
   0.400000000000000  34.396069793378309   2.072575255691516
   0.500000000000000  37.225819482896881   2.117788767225540
   0.600000000000000  40.267311806832083   2.176748354848179
   0.700000000000000  43.501170756497132   2.253434244955612
   0.800000000000000  46.935996455891825   2.350284334323296
   0.900000000000000  50.607214111441841   2.468336939929975
   1.000000000000000  54.530132303795753   2.610565196434151
   1.100000000000000  58.699942987824969   2.781877056174068
   1.200000000000000  63.091721492623776   2.989115289167727
   1.300000000000000  67.660426521509294   3.241057483112876
   1.400000000000000  72.340900152021618   3.548416043387018
   1.500000000000000  77.047867835923512   3.923838193047399
   1.600000000000000  81.675938399200817   4.381905972831026
   1.700000000000000  86.099604042062211   4.939136241154664
   1.800000000000000  90.173240338938967   5.613980674114787
   1.900000000000000  93.731106238485580   6.426825765487674
   2.000000000000000  96.587344063579224   7.399992826729324
   2.100000000000000  98.535979511319994   8.557737986975489
   2.200000000000000  99.305464073021838   9.923384989608529
   2.300000000000000  98.614325538390403  11.508508167056767
   2.400000000000000  96.285088095021379  13.306709884287148
   2.500000000000000  92.247162361659704  15.288207079976694
   2.600000000000000  86.585268003842529  17.394730789385640
   2.700000000000000  79.534876888241484  19.542682224883006
   2.800000000000000  71.536355306602971  21.622546215392298
   2.900000000000000  63.084804160958107  23.530018889724300
   3.000000000000000  54.623560359426598  25.181931216068016
   3.100000000000000  46.544096804423660  26.516262459375845
   3.200000000000000  39.186022392660291  27.492140181363510
   3.300000000000000  32.793221206846198  28.097843791494348
   3.400000000000000  27.336823355775415  28.376560995050070
   3.500000000000000  22.737483650943837  28.376423948829007
   3.600000000000000  18.913395131038023  28.142567798641700
   3.700000000000000  15.777073834688206  27.717805173699979
   3.800000000000000  13.235358800468109  27.142626186616940
   3.900000000000000  11.187311576411814  26.455624312391883
   4.000000000000000   9.527792494902240  25.691072088430193
   4.100000000000001   8.175847419742645  24.874012401250678
   4.200000000000000   7.068391701914113  24.024511571721689
   4.300000000000000   6.159197028246341  23.157992057962328
   4.400000000000000   5.418516351184799  22.285300010011660
   4.500000000000000   4.812949938211244  21.416178920107392
   4.600000000000001   4.312378013540997  20.558183397698095
   4.700000000000000   3.896728033888424  19.716433768708011
   4.800000000000001   3.549963138651240  18.894836312980246
   4.900000000000000   3.260082149910534  18.096083264276789
   5.000000000000000   3.018079494247696  17.321894302742177
   5.100000000000001   2.815390957478957  16.573545734604863
   5.200000000000000   2.645975745866385  15.851549329681150
   5.300000000000001   2.504682806034744  15.156183855662157
   5.400000000000000   2.387154825623396  14.487520111000972
   5.500000000000000   2.289828233286297  13.845420924912609
   5.600000000000001   2.209933198691999  13.229541157374021
   5.700000000000000   2.145493632523649  12.639327699124115
   5.800000000000001   2.095302528876663  12.074025260038512
   5.900000000000000   2.057819970371494  11.532952029583695
   6.000000000000000   2.031734456202639  11.015418711177421
   6.100000000000001   2.016223775353870  10.520656535937128
   6.200000000000000   2.010602361556739  10.047895063292874
   6.300000000000001   2.014321293290578   9.596362180987299
   6.400000000000000   2.026968293782502   9.165284105075676
   6.500000000000000   2.048267731007400   8.753885379925855
   6.600000000000001   2.078080617687949   8.361388878218298
   6.700000000000000   2.116388457892334   7.987017842929896
   6.800000000000001   2.163225001305634   7.630030231078177
   6.900000000000000   2.218783549199125   7.289725741835887
   7.000000000000000   2.283347851980680   6.965426060221426
   7.100000000000001   2.357279290720927   6.656474302505986
   7.200000000000000   2.441016877153240   6.362235016213561
   7.300000000000001   2.535077253673752   6.082094180120929
   7.400000000000000   2.640054693341341   5.815459204257671
   7.500000000000000   2.756621099877643   5.561758929906152
   7.600000000000000   2.885526007667043   5.320443629601537
   7.699999999999999   3.027596581756677   5.090985007131783
   7.800000000000000   3.183737617856438   4.872876197537634
   7.899999999999999   3.354931542338962   4.665631767112639
   8.000000000000000   3.542377769952171   4.468790794779848
   8.100000000000000   3.747750865739365   4.281929919898753
   8.199999999999999   3.972287358592226   4.104642734922401
   8.300000000000001   4.217388158593900   3.936541934626040
   8.399999999999999   4.484692559397574   3.777260394101559
   8.500000000000000   4.776078238226527   3.626451168757456
   8.600000000000000   5.093661255874061   3.483787494318877
   8.699999999999999   5.439796056703555   3.348962786827587
   8.800000000000001   5.817075468648450   3.221690642641979
   8.899999999999999   6.228330703212222   3.101704838437082
   9.000000000000000   6.676631355468448   2.988759331204541
   9.100000000000000   7.165285404060722   2.882628258252642
   9.199999999999999   7.697839211202703   2.783105937206295
   9.300000000000001   8.278077522678160   2.690006866007034
   9.399999999999999   8.910023467840816   2.603165722913034
   9.500000000000000   9.597978929632328   2.522479626174848
   9.600000000000000  10.346846934533836   2.447948500050266
   9.699999999999999  11.162039893459859   2.379385045866726
   9.800000000000001  12.049420139992524   2.316643469274441
   9.899999999999999  13.015307238019034   2.259644421438928
  10.000000000000000  14.066477981731850   2.208374999040992
  10.100000000000000  15.210166395628448   2.162888744276744
  10.199999999999999  16.454063734511507   2.123305644857589
  10.300000000000001  17.806318483488852   2.089812134010228
  10.399999999999999  19.275536357973380   2.062661090476664
  10.500000000000000  20.870780303683247   2.042171838514192
  10.600000000000000  22.601570496641628   2.028730147895408
  10.699999999999999  24.477884343176854   2.022788233908206
  10.800000000000001  26.510156479922557   2.024864757355776
  10.899999999999999  28.709278773817164   2.035544824556605
  11.000000000000000  31.084367983721659   2.055829025420224
  11.100000000000000  33.641600133885817   2.087482853795129
  11.199999999999999  36.398209815499314   2.130849364961110
  11.300000000000001  39.367967703938106   2.186873814289113
  11.400000000000000  42.559757627422172   2.257313701026955
  11.500000000000000  45.977576567015859   2.344738768299329
  11.600000000000000  49.620534656627697   2.452531003107799
  11.699999999999999  53.482855183010443   2.584884636330803
  11.800000000000001  57.553874585761008   2.746806142723651
  11.900000000000000  61.818042457320516   2.944114240918524
  12.000000000000000  66.254921542974230   3.183439893424474
  12.100000000000000  70.842835987138443   3.469191596026587
  12.199999999999999  75.522162277682469   3.812925452239469
  12.300000000000001  80.190042341363892   4.234758224058336
  12.400000000000000  84.721311371072218   4.754536876029922
  12.500000000000000  88.968565157515215   5.391807215404708
  12.600000000000000  92.762160089218654   6.165813892136889
  12.699999999999999  95.910213152526353   7.095500398884376
  12.800000000000001  98.198601931600223   8.199509071008816
  12.900000000000000  99.390964608420120   9.496181086575533
  13.000000000000000  99.228699962784120  11.003556466353615
  13.100000000000000  97.436046146590030  12.738369650117397
  13.199999999999999  93.897248393422203  14.679124234884609
  13.300000000000001  88.688453695875353  16.763830256803850
  13.400000000000000  82.012934272797054  18.915104969188761
  13.500000000000000  74.220377956306322  21.037723489086446
  13.600000000000000  65.797111775025812  23.020131785390436
  13.699999999999999  57.253976641417438  24.761724756250914
  13.800000000000001  49.023767724601868  26.194780217083739
  13.900000000000000  41.434244003664332  27.279068723563903
  14.000000000000000  34.706174448093314  28.001620983792861
  14.100000000000000  28.935101031065749  28.380190468466417
  14.199999999999999  24.065998570034740  28.461020682665570
  14.300000000000001  20.010527449455722  28.291941938864412
  14.400000000000000  16.675929778091504  27.918123986722730
  14.500000000000000  13.965024013772720  27.382080322488665
  14.600000000000000  11.774235192175343  26.723934672852753
  14.699999999999999   9.998177418711006  25.978837203372350
  14.800000000000001   8.554008385420849  25.173646553167053
  14.900000000000000   7.376179856212027  24.330003419116061
  15.000000000000000   6.415769580325549  23.464512074444599

火箭升空模型

火箭模型复现之后发现与课本数据有少许误差,我认为是matlab版本以及计算精度造成的差异

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
ts1 = 0:0.1:60;
ts2 = 60:0.1:71.5;
[t1,x1] = ode45(@(t,y)huojian1(t,y),ts1,[0,0]);%第一阶段
[t2,x2] = ode45(@(t,y)huojian2(t,y),ts2,[x1(length(x1),1),x1(length(x1),2)]);%第二阶段
figure(1);
plot(t1,x1(:,1),t2,x2(:,1));%速度图
%[t2,x2] = ode45(@huojian2,t2,[0 0]);
figure(2);
%画高度图
plot(t1,x1(:,2),t2,x2(:,2));
%画加速度图
a1 =diff(x1(:,1));
s1 = linspace(0,60,600)';%由于前向差分的缘故数据点与横坐标差了1,因此做一个数据拟合
s2 = linspace(60,71.5,115)';
a11 = spline(s1,a1,t1);
a2 = diff(x2(:,1));
a22 = spline(s2,a2,t2);
figure(3);
plot(t1,a11,t2,a22);grid on;
两阶段微分方程组
function dx = huojian1(t,x)
%阶段一
M = 1400;m = 18;F = 32000;k = 0.4;g = 9.8;
dx = [(F-k*x(1,1)*x(1,1))/(M-m*t)-g;x(1,1)];
end
function dx =huojian2(t,x)
%阶段二
k = 0.4;g = 9.8;M0=320;
dx = [(-k*x(1)*x(1))/M0-g;x(1)];
end

相关Blog
参考网址:https://blog.csdn.net/qq_29831163/article/details/89703911

$$