• 用MATLAB求解微分方程

    第一篇为 基础概念 ，第二篇为 R-K法的具体实现方法。

    （一）常微分方程的MATLAB求解

    概要：

    1. 常微分方程的MATLAB求解分为解析解、数值解
    2. 解析解(只有少数微分方程组有解析解)：dsolve函数
    3. 数值解：solver函数，MATLAB内部自带了许多求解器（各种数值求解的方法，如 R-K法）
    4. 常用函数文件构造微分方程组
    5. 数值求解时：高阶微分方程化成常微分方程组

    一、解析解

    dsolve函数：

    • y=dslove('eq1' , 'eq2' , ... , 'cond1' , 'cond2' , ... , 'v')
      • y——输出
      • eq1、eq2——微分方程
      • cond1、cond2——初值条件，省略时表示求通解
      • v——自变量，省略时默认为t
    • 微分方程中用 D 表示 自变量的导数
      • Dy—— y′
      • D2y—— y″
      • D3y—— y‴
    • 若找不到解析解，则返回其积分形式

    省略自变量和指定自变量的例子：

    1. >> dsolve('Dy = 2*x' , 'x')    %自变量x
    2. ans =
    3. x^2 + C1
    4.  
    5. >> dsolve('Dy = 2*x')    %默认自变量t
    6. ans =
    7. C2 + 2*t*x

    解析解例题：

    例1：求微分方程 dy/dx+2xy=xe^{-x^{2}}的通解

    1. >> y=dsolve('Dy+2*x*y = x*exp(-x^2)' , 'x')
    2. y =
    3. C3*exp(-x^2) + (x^2*exp(-x^2))/2

    例2：求下面微分方程的特解

    1. >> y=dsolve('D2y+4*Dy+29*y=0' , 'y(0)=0 ,Dy(0)=15' , 'x')
    2. y =
    3. 3*sin(5*x)*exp(-2*x)

    例3：求微分方程组的通解

    1. >> [x , y , z]=dsolve('Dx =2*x-3*y+3*z' , 'Dy=4*x-5*y+3*z' , 'Dz=4*x-4*y+2*z' , 't')
    2.  
    3. x =
    4. C6*exp(2*t) + C7*exp(-t)
    5.  
    6. y =
    7. C6*exp(2*t) + C7*exp(-t) + C8*exp(-2*t)
    8.  
    9. z =
    10. C6*exp(2*t) + C8*exp(-2*t)

    二、数值解

    只有很少一部分微分方程 (组) 能求出解析解；

    大部分微分方程（组）只能利用 数值方法 近似求解

    1、solver函数：

    • [T , Y] = solver(odefun , tspan , y0)
      • odefun——待求解的微分方程，可用matlab函数表示（m文件名）
      • tspan——求解区间
      • y0——初值条件(即自变量取最小值时（不一定为0）的函数值)
      • T——（向量）返回分割点的值，作自变量
      • Y——（向量）返回函数在分割点上的函数值
    • MATLAB在数值求解时 自动对求解区间进行分割
    • solver 为MATLAB的 ODE求解器（求解器中可以调用各种数值求解法：ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）
    求解器特点说明
    ode454，5阶 R-K 法首选方法
    ode232，3阶 R-K 法精度较低
    ode113Adams法，高低精度均可到 10^-3 ~10^-6计算时间比ode45短
    ode23t梯形算法（改进欧拉法）
    ode15s多步法：精度中等ode45失效时，可用
    ode23s单步法：低精度当精度较低时，计算时ode15s短
    ode23tb梯形算法：低精度当精度较低时，计算时比ode15s短

    ——————————————————————————————————

    2、odenfun（待求解的微分方程）的构建：

    （1）odenfun为 显示常微分方程（类似于显函数，可分离变量）

    • 可以用 inline函数 定义：直接将函数表达式内嵌在命令行里
      变量名=inline('函数表达式' , '变量名1' , '变量名2' , ...)
    • 也可以在 函数文件 中定义，然后通过函数句柄调用

    例1：求一阶常微分非线性方程的数值解，求解范围 [0,0.5]

    1. >> fun=inline('-2*y+2*x^2+2*x','x','y');	%定义常微分方程
    2. >> [x,y] = ode23(fun , [0,0.5] , 1)		%数值求解法使用ode23，自变量范围[0,0.5],步长默认,初值为1
    3.  
    4. x =
    5.          0
    6.     0.0400
    7.     0.0900
    8.     0.1400
    9.     0.1900
    10.     0.2400
    11.     0.2900
    12.     0.3400
    13.     0.3900
    14.     0.4400
    15.     0.4900
    16.     0.5000
    17.  
    18.  
    19. y =
    20.     1.0000
    21.     0.9247
    22.     0.8434
    23.     0.7754
    24.     0.7199
    25.     0.6764
    26.     0.6440
    27.     0.6222
    28.     0.6105
    29.     0.6084
    30.     0.6154
    31.     0.6179
    32.     
    33.     
    34.     
    35.     
    36. >> [x,y] = ode23(fun , [0:0.1:0.5] , 1)		%自定义步长为0.1
    37. x =
    38.          0
    39.     0.1000
    40.     0.2000
    41.     0.3000
    42.     0.4000
    43.     0.5000
    44.  
    45. y =
    46.     1.0000
    47.     0.8287
    48.     0.7103
    49.     0.6388
    50.     0.6093
    51.     0.6179

    ——————————————————————————————————————

    （2）求 高阶常微分方程，需将其化为 一阶常微分方程组 ，用函数文件定义

    函数文件格式 function dy = vdp1000(x , y)————其中x是向量，y是

    数值解例题：

    例：求解：

    解： 令 

    则原式子化成:

    vdp1000.m

    1. function dy = vdp1000(t , y)	%编写函数vdp1000,输入的t是自变量（不可省略），在后面的数值计算时，将自变量范围赋给t；y是一个2列的行向量，因为我们后面赋初值是有2个，且输出值也有2个
    2. dy = zeros(2,1);		%建立一个2行的列向量
    3. dy(1) = y(2);			%按照上边的式子赋值
    4. dy(2) = 1000*(1-y(1)^2)*y(2)-y(1);		%按照上边的式子赋值

    命令窗口

    1. >> [T , Y]=ode15s('vdp1000' , [0 , 3000] , [2 0]);		%调用方法，自变量t的范围[0 , 3000]，将初值 2,0 赋给 y(1),y(2) ,T——n个元素的列向量，Y——n行2列的矩阵，Y(:1)是式子中的y1，Y(:2)是式子中的y2
    2. >> plot(T , Y(:,1))		%显示Y(:1)，即式子中的y1，即原始式子的因变量x

    例2：手写 欧拉法 求常微分方程，x 范围 [0,0.5] , 步长 h=0.1

    解：根据欧拉公式： yn+1=yn+hf(xn,yn) -------------------(注：这里的 f(x,y)=y′ )

    代入 y′ , 可以得到 yn+1=0.9yn+0.1xn+0.1

    Euler.m

    1. function y=Euler()
    2.  
    3. x=[0:0.1:0.5];		%定义x
    4. y=zeros(size(x));	%创建y
    5. h=0.1;				%步长h=0.1
    6. len = length(x);	%向量x的长度
    7. y(1)=1;				%y初始化
    8. for i=2:len			%从y2开始计算
    9.     y(i)=0.9*y(i-1)+h*x(i-1) + 0.1;
    10. end
    11.  
    12. %画图
    13. figure
    14. hold on 
    15. grid on
    16. plot(x , y ,'r')		%红色线为数值解
    17. plot(x , x+exp(-x),'b')	%蓝色线为解析解
    18. xlabel('x')
    19. ylabel('y')
    20. title('Euler')

    命令窗口

    1. %函数调用
    2. >> y=Euler()
    3. y =
    4.     1.0000    1.0000    1.0100    1.0290    1.0561    1.0905

    例3：追击问题

    设位于坐标原点的甲舰向位于 x 轴上 A(1,0) 处的乙舰发射导弹，导弹头始终对准乙舰。

    已知：乙舰以最大的速度 v0 沿平行于 y 轴的直线行驶，导弹的速度是 5v0 。

    要求：模拟导弹的运行轨迹，并求击中点的坐标

    解：

    建立坐标系

    设 导弹的飞行曲线方程 为 y=y(x)

    在时刻 t 位于 P(x,y)

    eql.m

    1. function dy=eql(x,y)	%一阶微分方程组
    2. dy=zeros(2,1);		%构造2行的列向量
    3. dy(1)=y(2);				%dy(1)
    4. dy(2)=1/5*sqrt(1+y(2)^2)/(1-x);		dy(2)

    命令串口

    1. >> [x , y]=ode15s('eql',[0:0.001:1],[0 0]);		%数值求解微分方程组
    2. >> plot(x,y(:,1),'r.')					%画导弹轨迹，dy(1)，即原式子的y
    3. hold on
    4. >> y2=0:0.001:2;								
    5. >> plot(1,y2,'b*')					%画乙舰轨迹
    6. axis( [0 , 1 , 0 , 0.25 ] )				%坐标轴范围
    7. >> legend('导弹轨迹' , '乙舰轨迹')			%图例

    （二）常微分方程的MATLAB求解——手写龙格-库塔法（R-K法）模板、一般步骤 

    概要：

    • 龙格-库塔（R-K）法的写法：就是不断调用微分方程组，迭代计算出对于K1,K2，...，最后再叠加。需要注意的是高阶微分方程，其原函数的导数也是通过迭代计算得到的
    • 本人归纳了其套用 R-K 法的一般套路：3个函数、3个步骤——这也是MATLAB自带的求解方法的步骤
      • 三个函数：
        • Fun函数——用于存放一阶微分方程组
        • RK函数——用于使用几阶的R-K法求数值解，上边我只写了 4阶R-K法
        • 赋初值函数——只是单纯的写 x的范围，步长h，矩阵y的阶数，原函数的各个初值；以及调用 RK函数
      • 三个步骤：
        • 赋初值：写 x的范围，步长h，矩阵y的阶数，原函数的各个初值
        • 将高阶微分方程 拆分成 一阶微分方程组
        • 修改 Fun 函数：注意——向量dy的长度，和高阶微分方程的阶次有关
    • 本人也提供了 手写R-K法的模板，供大家使用参考

    例一：一阶微分方程

    使用MATLAB自带的 ode45函数 和 手写4阶R-K法 求解下列方程

    1.ode45函数解法：

    微分方程组

    1. function dy=Fun(x,y)
    2. dy=zeros(size(y));
    3.  
    4. dy(1) = sin(x)+y(1);					   %dy(1)表示以y的一阶导为f(x,y)

    调用ode45函数

    1. function y = MATLAB_RK()
    2. [X ,Y]=ode45('Fun' , [0 :0.1 :20] , [1]);		% y的初值1
    3.  
    4. %画图
    5. hold on
    6. grid on
    7. plot(X , Y(:,1))		%输出 Y的第一列，即原函数关于x的曲线
    8. xlabel('x')
    9. ylabel('y')

    2.手写4阶R-K解法：

    微分方程组

    1. function dy=Fun(x,y)
    2. dy=zeros(size(y));
    3.  
    4. dy(1) = sin(x)+y(1);					   %dy(1)表示以y的一阶导为f(x,y)

    调用手写R-K函数

    1. function y=PlotAll()
    2. % 准备阶段
    3. x=[0 : 0.1 :20];			% x——行向量
    4. h=0.1;					% 步长
    5. y=zeros(length(x) , 1);		        % y——矩阵，行数是x的元素个数，列数是几阶微分方程
    6. 					% 后面的1表示一阶微分方程,如果是2则表示二阶微分方程，以此类推
    7. y(1,:)=[1];      			% y，第一行赋初值，几阶微分方程就要赋几个初值
    8.  
    9. % 函数调用
    10. Y=FourRK(x,h,y);			% 调用4阶R-K法，输入参数x，h，y
    11.  
    12.  
    13. %画图
    14. hold on
    15. grid on
    16. plot(x , Y(:,1))		%输出 Y的第一列，即原函数关于x的曲线
    17. xlabel('x')
    18. ylabel('y')

    4阶R-K函数

    1. function y=FourRK(x,h,y)
    2.  
    3. len = length(x);	
    4.  
    5. for i=2:len		%循环：直到求完最后一个x取值
    6.  
    7. K1 = Fun(x(i-1),y(i-1,:));  		% K1
    8.  
    9. K2 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K1);		% K2
    10.  
    11. K3 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K2);		% K3
    12.  
    13. K4 = Fun(x(i-1)+h , y(i-1,:)+h.*K3);			% K4
    14.  
    15.  
    16. y(i,:) = y(i-1,:)+h/6.*(K1 + 2*K2 + 2*K3 + K4);           % 更新y的第i行，y(i,1)是原函数的y，y(i,2)是原函数的一阶导，以此类推
    17. 														  % 因为这里是1阶微分方程，所以只有y(i,1)，无y(i,2)
    18. end

    例二：二阶微分方程

    使用MATLAB自带的 ode45函数 和 手写4阶R-K法 求解下列方程

    解：将高阶微分方程 化成 一阶微分方程组

    我们下面的 Fun函数 就是依据这个一阶微分方程组写的

    1.ode45函数解法：

    微分方程组

    1. function dy=Fun(x,y)			  % y(1)表示原函数y，y(2)表示原函数y的一阶导
    2. dy=zeros(size(y));
    3. dy(1)=y(2);                               % dy(1)表示以y的一阶导为f(x,y)
    4. dy(2)=-1/2*(y(1)+(y(1)*y(1)-1)*y(2));     % dy(2)表示以y的二阶导为f(x,y)

    调用ode45函数

    1. function y = MATLAB_RK()
    2. [X ,Y]=ode45('Fun' , [0 :0.1 :20] , [0.25 , 0.0]);		% y的初值0.25 ，y一阶导的初值 0.0
    3.  
    4. %画图
    5. hold on
    6. grid on
    7. plot(X , Y(:,1))		%输出 Y的第一列，即原函数关于x的曲线
    8. xlabel('x')
    9. ylabel('y')

    2.手写4阶R-K解法：

    微分方程组

    1. function dy=Fun(x,y)
    2. dy=zeros(size(y));
    3. dy(1)=y(2);                               %dy(1)表示以y的一阶导为f(x,y)
    4. dy(2)=-1/2*(y(1)+(y(1)*y(1)-1)*y(2));     %dy(2)表示以y的二阶导为f(x,y)

    调用手写R-K函数

    1. function y=PlotAll()
    2. % 准备阶段
    3. x=[0 : 0.1 :20];		% x——行向量
    4. h=0.1;				% 步长
    5. y=zeros(length(x) , 2);	        % y——矩阵，行数是x的元素个数，列数是几阶微分方程
    6. 				% 后面的1表示一阶微分方程,如果是2则表示二阶微分方程，以此类推。这里是二阶微分方程，所以是2
    7. y(1,:)=[0.25 , 0.0];            % y，第一行赋初值，几阶微分方程就要赋几个初值。
    8. 				% y的初值0.25 ，y一阶导的初值 0.0
    9.  
    10. % 函数调用
    11. Y=FourRK(x,h,y);		% 调用4阶R-K法，输入参数x，h，y
    12.  
    13. %画图
    14. hold on
    15. grid on
    16. plot(x , Y(:,1))		%输出 Y的第一列，即原函数关于x的曲线
    17. xlabel('x')
    18. ylabel('y')

    4阶R-K函数

    1. function y=FourRK(x,h,y)
    2.  
    3. len = length(x);
    4.  
    5. for i=2:len			%循环：直到求完最后一个x取值
    6.  
    7. K1 = Fun(x(i-1),y(i-1,:));  		% K1
    8.  
    9.  
    10. K2 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K1);		% K2
    11.  
    12.  
    13. K3 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K2);		% K3
    14.  
    15.  
    16. K4 = Fun(x(i-1)+h , y(i-1,:)+h.*K3);		% K4
    17.  
    18.  
    19. y(i,:) = y(i-1,:)+h/6.*(K1 + 2*K2 + 2*K3 + K4);           % 更新y的第i行，y(i,1)是原函数的y，y(i,2)是原函数的一阶导，以此类推
    20. 														  % 因为这里是2阶微分方程，所以有y(i,1)，有y(i,2)
    21. end
    22.  
    23.  
    24.  
    25.  
    26.  
    27. %本质上，将高阶微分化为一阶微分方程组，然后按照各自当作f(x,y），分别进行迭代，最后更新

    例二 手写R-K法 答疑：

    一阶微分方程可能还好，二阶微分方程可能就没那么好理解了，为此对于例二我换了种 不通用但更直观 的写法，帮助大家理解

    首先得理解为什么要将 二阶微分方程 拆成 一阶微分方程组 ，因为这样我们就可以根据 设 f(x,y)=y 的一阶导数，进行迭代求解。————如此题，我们拆成了2个一阶微分方程，所以就需要就同时进行2个R-K迭代

    • 这里我将 y 和 y' 拆开成2个行向量——————上边是 将2个合并到一个矩阵 y 中， y(1)为原函数y，y(2)为原函数y的一阶导
    • 并且我们需要分别定义函数，来迭代计算 y 和 y' ——————上边是 将2个合并到一个矩阵 y 中.一起计算了

    Fun1————迭代求 y

    1. function f1=Fun1(x,y,y1)
    2.  
    3. f1=y1;      %f2表示以y的一阶导为f(x,y)——f(x,y)=输入的y1

    Fun2————迭代求 y‘

    1. function f2=Fun2(x,y,y1)
    2.  
    3. f2=-1/2*(y+(y*y-1)*y1);     %f2表示以y的二阶导为f(x,y)

    4阶R-K函数

    1. function y=FourRK2()
    2.  
    3. x=[0 : 0.1 :20];
    4. y=zeros(size(x));
    5. y1=zeros(size(x));
    6. h=0.1;
    7. len = length(x);
    8. y(1)=0.25;      % 原函数y初值
    9. y1(1)=0.0;      % 原函数y的一阶导初值
    10.  
    11. for i=2:len
    12. %Kmn 表示第n次迭代，以m次导数作为f(x,y)
    13.  
    14. K11 = Fun1(x(i-1),y(i-1),y1(i-1));  
    15. K21 = Fun2(x(i-1),y(i-1),y1(i-1));
    16.  
    17. K12 = Fun1(x(i-1)+1/2*h , y(i-1)+1/2*h*K11 , y1(i-1)+1/2*h*K21);
    18. K22 = Fun2(x(i-1)+1/2*h , y(i-1)+1/2*h*K11 , y1(i-1)+1/2*h*K21);
    19.  
    20. K13 = Fun1(x(i-1)+1/2*h , y(i-1)+1/2*h*K12 , y1(i-1)+1/2*h*K22);
    21. K23 = Fun2(x(i-1)+1/2*h , y(i-1)+1/2*h*K12 , y1(i-1)+1/2*h*K22);
    22.  
    23. K14 = Fun1(x(i-1)+h , y(i-1)+h*K13 , y1(i-1)+h*K23);
    24. K24 = Fun2(x(i-1)+h , y(i-1)+h*K13 , y1(i-1)+h*K23);
    25.  
    26. y(i) = y(i-1)+h/6*(K11 + 2*K12 + 2*K13 + K14);           % 更新下一个原函数y
    27. y1(i) = y1(i-1)+h/6*(K21 + 2*K22 + 2*K23 + K24);         % 更新下一个原函数y的一阶导
    28. end
    29.  
    30.  
    31. %画图
    32. hold on
    33. grid on
    34. plot(x , y)		
    35. xlabel('x')
    36. ylabel('y')
    37.  
    38.  
    39. %本质上，将高阶微分化为一阶微分方程组，然后按照各自当作f(x,y），分别进行迭代，最后更新

    总结：手写微分方程的 模板 和 一般步骤

    一、原理讲解：

    1.模板：

    • Fun函数——用于存放一阶微分方程组
    • RK函数——用于使用几阶的R-K法求数值解，上边我只写了 4阶R-K法
    • 赋初值函数——只是单纯的写 x的范围，步长h，矩阵y的阶数，原函数的各个初值；以及调用 RK函数

    2.一般步骤：

    其实这里 手写的R-K法 原理和步骤就是 MATLAB自带的求解R-K方法

    • 赋初值：写 x的范围，步长h，矩阵y的阶数，原函数的各个初值
    • 将高阶微分方程 拆分成 一阶微分方程组
    • 修改 Fun 函数：注意——向量dy的长度，和高阶微分方程的阶次有关

    二、具体模板：

    1.微分方程模板和准备函数模板

    Fun函数————需要自己修改

    1. function dy=Fun(x,y)
    2. dy=zeros(size(y));
    3. dy(1)=    ;                               %dy(1)表示以y的一阶导为f(x,y)
    4. dy(2)=    ;     			  %dy(2)表示以y的二阶导为f(x,y)
    5. dy(3)=    ;     			  %dy(3)表示以y的三阶导为f(x,y)
    6. ...

    准备初值函数————需要自己修改值，需要修改的用【】括起了

    1. function Y=PlotAll()
    2. % 准备阶段
    3. x=[0 : 0.1 :20];		% 【x——行向量】
    4. h=0.1;				% 【步长】
    5. y=zeros(length(x) , 2);	        % 【y——矩阵，行数是x的元素个数，列数是几阶微分方程】
    6. 				% 后面的1表示一阶微分方程,如果是2则表示二阶微分方程，以此类推。这里是二阶微分方程，所以是2
    7. y(1,:)=[0.25 , 0.0];            % 【y，第一行赋初值，几阶微分方程就要赋几个初值。】
    8. 				% 【y的初值0.25 ，y一阶导的初值 0.0】
    9.  
    10. % 函数调用
    11. Y=FourRK(x,h,y);		% 调用4阶R-K法，输入参数x，h，y
    12.  
    13. %画图
    14. hold on
    15. grid on
    16. plot(x , Y(:,1))		%输出 Y的第一列，即原函数关于x的曲线
    17. xlabel('x')
    18. ylabel('y')

    2.欧拉法、改进欧拉法、4阶R-K法模板

    1阶R-K函数(欧拉法)

    1. function y=OneRK(x,h,y)
    2.  
    3. len = length(x);
    4.  
    5. for i=2:len			%循环：直到求完最后一个x取值
    6.  
    7. K1 = Fun(x(i-1),y(i-1,:));
    8.  
    9. y(i,:)=y(i-1,:)+h.*K1;		% 更新y的第i行，y(i,1)是原函数的y，y(i,2)是原函数的一阶导，以此类推
    10. 	
    11. end
    12.  
    13.  
    14.  
    15.  
    16.  
    17. %本质上，将高阶微分化为一阶微分方程组，然后按照各自当作f(x,y），分别进行迭代，最后更新

    2阶R-K函数(改进欧拉法)

    1. function y=TwoRK(x,h,y)
    2.  
    3. len = length(x);
    4.  
    5. for i=2:len			%循环：直到求完最后一个x取值
    6.  
    7. K1 = Fun(x(i-1) , y(i-1,:));
    8. K2 = Fun(x(i-1)+h , y(i-1,:)+h.*K1);
    9.  
    10. y(i,:)=y(i-1,:)+h/2.*(K1+K2);		% 更新y的第i行，y(i,1)是原函数的y，y(i,2)是原函数的一阶导，以此类推
    11. end
    12.  
    13.  
    14.  
    15.  
    16.  
    17. %本质上，将高阶微分化为一阶微分方程组，然后按照各自当作f(x,y），分别进行迭代，最后更新

    4阶R-K函数

    1. function y=FourRK(x,h,y)
    2.  
    3. len = length(x);
    4.  
    5. for i=2:len			%循环：直到求完最后一个x取值
    6.  
    7. K1 = Fun(x(i-1),y(i-1,:));  		% K1
    8.  
    9.  
    10. K2 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K1);		% K2
    11.  
    12.  
    13. K3 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K2);		% K3
    14.  
    15.  
    16. K4 = Fun(x(i-1)+h , y(i-1,:)+h.*K3);		% K4
    17.  
    18.  
    19. y(i,:) = y(i-1,:)+h/6.*(K1 + 2*K2 + 2*K3 + K4);           % 更新y的第i行，y(i,1)是原函数的y，y(i,2)是原函数的一阶导，以此类推
    20. end
    21.  
    22.  
    23.  
    24. %本质上，将高阶微分化为一阶微分方程组，然后按照各自当作f(x,y），分别进行迭代，最后更新

    补充：关于解为lambertw的微分方程数值解

    例，求解下面的微分方程

    可以发现这个初值会让 y′(0) 的值有无数个。如果我们写 dy=xy/(1−y) 就会出现分母为0的情况，一开始我就是这么定义 Fun 函数的，发现结果为 NaN。

    正解：

    • 其实这个函数与 y′ 没有关系，不管 y′(0) 取多少，函数图像总是先为0，后面快速趋向正无穷
    • 因此我们就默认取 dy=0 (当然大家也可以取其他值)，然后不将 dy 到一边，直接用原表达式计算

    Fun

    1. function dy=Fun(x,y)
    2.  
    3. dy=-200;
    4. dy=(x*y)+dy*y;

  • 相关阅读:
    高压配电柜电力安全监控解决方案
    uni-app开发android应用流程
    java计算机毕业设计志愿者社会服务管理系统源码+mysql数据库+系统+lw文档+部署
    【1805. 字符串中不同整数的数目】
    探索MySQL错误: 1241 - Operand should contain 1 column(s)问题解决方案
    JS本地存储技术
    智能优化算法——混合领导优化算法（Matlab&Matlab代码实现）
    拯救SQL Server数据库事务日志文件损坏的终极大招
    微信小程序：单行输入和多行输入组件
    推文项目进展如何
  • 原文地址：https://blog.csdn.net/sinat_41942180/article/details/136535223