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

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

$\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)$

例子:
$\{\begin{array}{l}\frac{dy}{dx}=xy\\ y_0=1\end{array}$

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

如果通过一阶差分求解
$\{\begin{array}{l}{y}_{n+1}=y_n+h{x}_n{y}_n\\ x_{n+1}=x_n+h\\ y_0=1\\ x_0=0\end{array}$

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]);

如果再加上解析解

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]);

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]);

细看是两条线

应用:小船渡河问题:

河宽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}$

$\{\begin{array}{l}\frac{dy}{dt}=V_{2}si{n}_{\theta }\\ \frac{dx}{dt}=X_1+V_{2}co{s}_{\theta }\end{array}$

$\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]);

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

进阶求二阶微分方程

例子:
$\{\begin{array}{l}{y}^{\prime \prime }-2y^{\prime }+y=0\\ y_0=1\\ y_0^{\prime }=2\end{array}$
解析解为:$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]);