【数学建模】微分方程的数值求解

一阶差分求解微分方程原理:

$\frac{dy}{dx}=f(x_n,y_n)$

$\frac{y_{n+1}-y_n}{x_{n+1}-x_n}=f(x_n,y_n)$

令$h=x_{n+1}-x_n$

$y_{n+1}=y_n+hf(x_n,y_n)$

例子:
{ d y d x = x y y 0 = 1

$$\begin{cases} \dfrac{dy}{dx}=xy\\ y_0 = 1 \end{cases}$$ ⎩ ⎨ ⎧​dxdy​=xyy0​=1​

通过高等数学可以求得解析解为:$y=e^{\frac{1}{2}{x}^2}$

如果通过一阶差分求解
{ y n + 1 = y n + h x n y n x n + 1 = x n + h y 0 = 1 x 0 = 0

$$\begin{cases} y_{n+1}=y_n+hx_ny_n\\ x_{n+1}=x_n+h\\ y_0 = 1\\ x_0 = 0 \end{cases}$$ ⎩ ⎨ ⎧​yn+1​=yn​+hxn​yn​xn+1​=xn​+hy0​=1x0​=0​

matlab求解:

y = zeros(1,1000);
x = zeros(1,1000);
y(1)=1;
x(1)=0;
h=0.1;

for n=1:1000
    y(n+1)=y(n)+h*x(n)*y(n);
    x(n+1)=x(n)+h;
end
plot(x,y,'-');
xlim([0,6]);
ylim([0,100]);
1
2
3
4
5
6
7
8
9
10
11
12
13

如果再加上解析解

y = zeros(1,1000);
x = zeros(1,1000);
y(1)=1;
x(1)=0;
h=0.1;

for n=1:1000
    y(n+1)=y(n)+h*x(n)*y(n);
    x(n+1)=x(n)+h;
end
plot(x,y,'-');
hold on;
ezplot('exp(x^2/2)');
xlim([0,6]);
ylim([0,100]);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

h=0.01如下图

四阶龙格-库塔方法

原理/公式

优点:精度高

ezplot('exp(x^2/2)');
hold on;
%四阶龙格-库塔方法
y = zeros(1,1000);
x = zeros(1,1000);
y(1)=1;
x(1)=0;
h=0.1;
for n=1:1000
    k1 = h*(x(n)*y(n));%h*f(x(n),y(n))
    k2 = h*((x(n)+h/2)*(y(n)+k1/2));%h*f(x(n)+h/2,y(n)+k1/2)
    k3 = h*((x(n)+h/2)*(y(n)+k2/2));%h*f(x(n)+h/2,y(n)+k2/2)
    k4 = h*((x(n)+h)*(y(n)+k3));%h*f(x(n)+h,y(n)+k3)
    y(n+1)=y(n)+(1/6)*(k1+2*k2+2*k3+k4);
    y(1)=1;
    x(n+1)=x(n)+h;
end
plot(x,y,'-');
xlim([0,6]);
ylim([0,100]);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

细看是两条线

应用:小船渡河问题:

河宽100米，A(0,100) , O(0,0),水流速$V_1$,小船速度$V_2$,水流速度始终垂直OA并向下流冲（方向不变）
从O到A点，求轨迹路线图和时间
与$V_2$平行的相量$(-x,100,-y)$

$co{s}_{\theta }=\frac{-x}{\sqrt{x^2}+(100-y{)}^2}$

$si{n}_{\theta }=\frac{100-y}{\sqrt{x^2}+(100-y{)}^2}$

{ d y d t = V 2 s i n θ d x d t = X 1 + V 2 c o s θ

$$\begin{cases} \dfrac{dy}{dt}=V_2sin_{\theta}\\ \dfrac{dx}{dt}=X_1+V_2cos_{\theta} \end{cases}$$ ⎩ ⎨ ⎧​dtdy​=V2​sinθ​dtdx​=X1​+V2​cosθ​​

$\frac{dy}{dx}=\frac{V_{2}si{n}_{\theta }}{V_1+V_{2}co{s}_{\theta }}=\frac{V_2(100-y)}{\sqrt{x^2+(100-y{)}^2}{V}_1-V_{2}x}$

因为$x$变量随着$y$增大有增有减，而$y$变量随着$x$增大有递增，所以$y$为自变量可以求解

变换后:
$\frac{dx}{dy}=\frac{V_1+V_{2}co{s}_{\theta }}{V_{2}si{n}_{\theta }}=\frac{\sqrt{x^2+(100-y{)}^2}{V}_1-V_{2}x}{V_2(100-y)}$

$x_0=y_0=0$

y = zeros(1,1000);
x = zeros(1,1000);
x(1)=0;
y(1)=0;
v1=30;
v2=200;
h=0.1;
%差分求解轨迹
for n=1:1000 %1000*h = 100
    x(n+1)=x(n)+h*((sqrt(x(n)^2+(100-y(n))^2)*v1-v2*x(n))/(v2*(100-y(n))));
    y(n+1)=y(n)+h;
end

plot(x,y,'-');
hold on;
%四阶龙格-库塔方法求轨迹
y = zeros(1,1000);
x = zeros(1,1000);
x(1)=0;
y(1)=0;
for n=1:1000
    k1 = h*((sqrt(x(n)^2+(100-y(n))^2)*v1-v2*x(n))/(v2*(100-y(n))));%h*f(x(n),y(n))
    k2 = h*((sqrt((x(n)+h/2)^2+(100-(y(n)+k1/2))^2)*v1-v2*(x(n)+h/2))/(v2*(100-(y(n)+k1/2))));%h*f(x(n)+h/2,y(n)+k1/2)
    k3 = h*((sqrt((x(n)+h/2)^2+(100-(y(n)+k2/2))^2)*v1-v2*(x(n)+h/2))/(v2*(100-(y(n)+k2/2))));%h*f(x(n)+h/2,y(n)+k2/2)
    k4 = h*((sqrt((x(n)+h)^2+(100-(y(n)+k3))^2)*v1-v2*(x(n)+h))/(v2*(100-(y(n)+k3))));%h*f(x(n)+h,y(n)+k3)
    x(n+1)=x(n)+(1/6)*(k1+2*k2+2*k3+k4);
    y(n+1)=y(n)+h;
end
plot(x,y,'-');
%(sqrt(x(n)^2+(100-y(n))^2)
%差分求解时间
t = zeros(1,1000);
for n=1:1000
    t(n+1)=t(n)+h*(1/(v2*(100-y(n))/(sqrt(x(n)^2+(100-y(n))^2))));
end

xlim([0,10]);
ylim([0,100]);
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

使用特解验证时间算的对不对

进阶求二阶微分方程

例子:
{ y ′ ′ − 2 y ′ + y = 0 y 0 = 1 y 0 ′ = 2

$$\begin{cases} y''-2y'+y=0\\ y_0=1\\ y'_0=2 \end{cases}$$ ⎩ ⎨ ⎧​y′′−2y′+y=0y0​=1y0′​=2​
解析解为:$y=e^x+x{e}^x$
差分求解:
设$y_1=y,y_2=y^{\prime }$
则
$\{\begin{array}{l}{y}_1^{\prime }=y2\\ y_2^{\prime }=2y_2-y_1\\ y_{1_0}=1\\ y_{2_0}=2\end{array}$

y1 = zeros(1,1000);
y2 = zeros(1,1000);
x = zeros(1,1000);
y1(1)=1;
y2(1)=2;
x(1)=0;
h=0.01;
xlim([0,6]);
ylim([0,100]);

for n=1:1000
    y2(n+1)=y2(n)+h*(2*y2(n)-y1(n));
    y1(n+1)=y1(n)+h*y2(n);
    x(n+1)=x(n)+h;
end

%plot(x,y1,'-');
hold on;
%四阶龙格-库塔方法
y1 = zeros(1,1000);
y2 = zeros(1,1000);
x = zeros(1,1000);
y1(1)=1;
y2(1)=2;
x(1)=0;
for n=1:1000
    k1 = h*(2*y2(n)-y1(n));%h*f(x(n),y(n))
    k2 = h*(2*(y2(n)+k1/2)-(y1(n)+h/2));%h*f(x(n)+h/2,y(n)+k1/2)
    k3 = h*(2*(y2(n)+k2/2)-(y1(n)+h/2));%h*f(x(n)+h/2,y(n)+k2/2)
    k4 = h*(2*(y2(n)+k3)-(y1(n)+h));%h*f(x(n)+h,y(n)+k3)

    y2(n+1)=y2(n)+(1/6)*(k1+2*k2+2*k3+k4);
    
    k1 = h*y2(n);%h*f(x(n),y1(n))
    k2 = h*(y2(n)+h/2);%h*f(x(n)+h/2,y(n)+k1/2)
    k3 = h*(y2(n)+h/2);%h*f(x(n)+h/2,y(n)+k2/2)
    k4 = h*(y2(n)+h);%h*f(x(n)+h,y(n)+k3)
    y1(n+1)=y1(n)+(1/6)*(k1+2*k2+2*k3+k4);
    
    x(n+1)=x(n)+h;
end
plot(x,y1,'-');

%解析解

ezplot('exp(x)+x*exp(x)');
xlim([0,6]);
ylim([0,100]);
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