数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程

matlab解二阶微分方程组在数学和工程领域中，二阶微分方程组是一类常见的问题，它们描述了许多自然现象和工程系统的动力学行为。

在本文中，我们将讨论如何使用Matlab解二阶微分方程组，并探讨其中的一些技巧和注意事项。

让我们来看一个简单的例子。

假设我们有一个二阶微分方程组，描述了一个弹簧-质量-阻尼系统的运动。

我们可以使用Matlab的ode45函数来求解这个方程组，该函数是一种常用的数值求解器，可以有效地解决微分方程组的数值解。

在使用ode45函数时，我们需要定义一个包含微分方程的函数，并将其作为输入传递给ode45函数。

在定义这个函数时，我们需要注意输入参数的顺序和返回值的格式，以确保程序能够正确运行。

我们还需要设置初始条件和求解的时间范围。

通过调用ode45函数并传递这些参数，我们可以得到微分方程组的数值解，并将其存储在一个数组中以供进一步分析和可视化。

当然，对于复杂的二阶微分方程组，可能需要更多的技巧和方法来求解。

在这种情况下，我们可以考虑使用符号计算工具箱来分析微分方程的性质和解析解。

然后，我们可以将这些解析解转化为数值解，并进行比较和验证。

我们还可以利用Matlab的图形化界面工具来直观地展示微分方程组的解，并进行参数的调整和实时可视化。

这种方法可以帮助我们更好地理解系统的动态行为，并优化系统的设计和控制策略。

总的来说，使用Matlab解二阶微分方程组是一种强大而灵活的工具，在数学建模和工程应用中发挥着重要作用。

通过合理地选择数值求解方法、调整参数和优化算法，我们可以高效地求解复杂的微分方程组，并获得准确的数值解。

希望本文能够帮助读者更好地理解和应用Matlab在二阶微分方程组求解中的技巧和方法。

MATLAB中龙格库塔（RUNGEKUTTA）方法原理及实现函数功能ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kutta算法。

ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(Δx)3。

解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解.使用方法[T,Y]=ode45(odefun,tspan,y0)odefun是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名tspan是区间[t0tf]或者一系列散点[t0,t1,...,tf]y0是初始值向量T返回列向量的时间点Y返回对应T的求解列向量[T,Y]=ode45(odefun,tspan,y0,options)options是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等[T,Y,TE,YE,IE]=ode45(odefun,tspan,y0,options)在设置了事件参数后的对应输出TE事件发生时间YE事件解决时间IE事件消失时间sol=ode45(odefun,[t0tf],y0...)sol结构体输出结果应用举例1求解一阶常微分方程程序:一阶常微分方程odefun=@(t,y)(y+3*t)/t^2;%定义函数tspan=[14];%求解区间y0=-2;%初值[t,y]=ode45(odefun,tspan,y0);plot(t,y)%作图title('t^2y''=y+3t,y(1)=-2,1<t<4')< p="">legend('t^2y''=y+3t')xlabel('t')ylabel('y')%精确解%dsolve('t^2*Dy=y+3*t','y(1)=-2')%ans=一阶求解结果图%(3*Ei(1)-2*exp(1))/exp(1/t)-(3*Ei(1/t))/exp(1/t)2求解高阶常微分方程关键是将高阶转为一阶,odefun的书写.F(y,y',y''...y(n-1),t)=0用变量替换,y1=y,y2=y'...注意odefun方程定义为列向量dxdy=[y(1),y(2)....]程序:function Testode45tspan=[3.94.0];%求解区间y0=[28];%初值[t,x]=ode45(@odefun,tspan,y0);plot(t,x(:,1),'-o',t,x(:,2),'-*')legend('y1','y2')title('y''''=-t*y+e^t*y''+3sin2t')xlabel('t')ylabel('y')function y=odefun(t,x)y=zeros(2,1);%列向量y(1)=x(2);y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);endend高阶求解结果图相关函数ode23,ode45,ode113,ode15s,ode23s,ode23t,ode23tbMatlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。

标题：探究matlab中的变系数微分方程及其数值解法——龙格库塔在matlab中，我们常常会遇到变系数微分方程及其数值解法的问题。

变系数微分方程是一类常见的微分方程，在实际问题中具有重要的应用价值。

而龙格库塔是一种常用的数值解法，能够有效地求解变系数微分方程，具有广泛的适用性。

1. 变系数微分方程的概念与应用变系数微分方程是指微分方程中的系数是随自变量或未知函数的变量的而变化的微分方程，通常可以表示为形式：\( \frac{d}{dx}(a(x)y') + b(x)y = f(x) \)。

这类微分方程在物理、工程和经济学等领域中具有广泛的应用，例如在电路中的电流变化、非线性振动系统的描述以及经济增长模型等方面都有相关的应用。

2. 龙格库塔法的原理和适用性龙格库塔法是一种常用的数值解法，用于求解常微分方程组及其初值问题。

与欧拉法相比，龙格库塔法能够更准确地求解微分方程，并具有更好的数值稳定性和收敛性。

具体而言，龙格库塔法通过多步逼近的方式，利用对目标函数在不同点的取值进行加权平均，得到更精确的数值解。

3. 在matlab中求解变系数微分方程的实现在matlab中，我们可以利用ode45、ode23等内置函数，或者自行编写龙格库塔法的求解程序，来求解变系数微分方程。

其中，ode45函数是matlab中最常用的求解微分方程的函数之一，它基于龙格库塔法，能够高效地求解一阶和二阶的常微分方程，并具有良好的数值稳定性和收敛性。

4. 个人观点和总结在我看来，掌握变系数微分方程及其数值解法对于理解和应用微分方程具有重要意义。

通过学习龙格库塔法，我们不仅能够更深入地理解微分方程的数值解法，还能够应用于实际问题中，解决复杂的变系数微分方程。

我建议在学习matlab时，要深入了解变系数微分方程及其数值解法，特别是龙格库塔法，这将对未来的科研和工程实践大有裨益。

在本文中，我们探讨了matlab中的变系数微分方程及其数值解法，重点介绍了龙格库塔法在求解变系数微分方程中的应用。

MATLAB常微分⽅程数值解——欧拉法、改进的欧拉法与四阶龙格库塔⽅法MATLAB常微分⽅程数值解作者：凯鲁嘎吉 - 博客园1.⼀阶常微分⽅程初值问题2.欧拉法3.改进的欧拉法4.四阶龙格库塔⽅法5.例题⽤欧拉法，改进的欧拉法及4阶经典Runge-Kutta⽅法在不同步长下计算初值问题。

步长分别为0.2,0.4,1.0.matlab程序：function z=f(x,y)z=-y*(1+x*y);function R_K(h)%欧拉法y=1;fprintf('欧拉法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K=f(x,y);y=y+h*K;fprintf('欧拉法：x=%f, y=%f\n',x+h,y);endfprintf('\n');%改进的欧拉法y=1;fprintf('改进的欧拉法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h,y+h*K1);y=y+(h/2)*(K1+K2);fprintf('改进的欧拉法：x=%f, y=%f\n',x+h,y);endfprintf('\n');%龙格库塔⽅法y=1;fprintf('龙格库塔法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h/2,y+(h/2)*K1);K3=f(x+h/2,y+(h/2)*K2);K4=f(x+h,y+h*K3);y=y+(h/6)*(K1+2*K2+2*K3+K4);fprintf('龙格库塔法：x=%f, y=%f\n',x+h,y);end结果：>> R_K(0.2)欧拉法：x=0.000000, y=1.000000欧拉法：x=0.200000, y=0.800000欧拉法：x=0.400000, y=0.614400欧拉法：x=0.600000, y=0.461321欧拉法：x=0.800000, y=0.343519欧拉法：x=1.000000, y=0.255934改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=0.200000, y=0.807200改进的欧拉法：x=0.400000, y=0.636118改进的欧拉法：x=0.600000, y=0.495044改进的欧拉法：x=0.800000, y=0.383419改进的欧拉法：x=1.000000, y=0.296974龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=0.200000, y=0.804636龙格库塔法：x=0.400000, y=0.631465龙格库塔法：x=0.600000, y=0.489198龙格库塔法：x=0.800000, y=0.377225龙格库塔法：x=1.000000, y=0.291009>> R_K(0.4)欧拉法：x=0.000000, y=1.000000欧拉法：x=0.400000, y=0.600000欧拉法：x=0.800000, y=0.302400改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=0.400000, y=0.651200改进的欧拉法：x=0.800000, y=0.405782龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=0.400000, y=0.631625龙格库塔法：x=0.800000, y=0.377556>> R_K(1)欧拉法：x=0.000000, y=1.000000欧拉法：x=1.000000, y=0.000000改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=1.000000, y=0.500000龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=1.000000, y=0.303395注意：在步长h为0.4时，要将for i=1:1/h改为for i=1:0.8/h。

高树磊 2010010909题2--3%欧拉法0.1步长下输出响应曲线和真实响应曲线的比较subplot(221)x=1;%步长为0.1y=x;h=0.1;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'r')hold onx=1;y=x;h=0.1;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'*g')title('欧拉法0.1步长下输出响应曲线和真实响应曲线的比较'); xlabel('t/time');ylabel('y/gain');legend('Euler 0.1步长''真实输出');%欧拉法0.05步长下输出响应曲线和真实响应曲线的比较subplot(222)x=1;%步长为0.05y=x;h=0.05;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'k')hold onx=1;y=x;h=0.1;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'*y')title('欧拉法0.05步长下输出响应曲线和真实响应曲线的比较'); xlabel('t/time');ylabel('y/gain');legend('Euler 0.05步长','真实输出');%欧拉法0.5步长下输出响应曲线和真实响应曲线的比较subplot(223)x=1;%步长为0.5y=x;h=0.5;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'b')hold onx=1;y=x;h=0.1;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'*r')title('欧拉法0.5步长下输出响应曲线和真实响应曲线的比较'); xlabel('t/time');ylabel('y/gain');legend('Euler 0.5步长','真实输出');%欧拉法不同步长下输出响应曲线和真实响应曲线的比较subplot(224)x=1;%步长为0.1y=x;h=0.1;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'k')hold onx=1;%步长为0.05y=x;h=0.05;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'y')hold onx=1;%步长为0.5y=x;h=0.5;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'b')hold onx=1;y=x;h=0.1;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'*r')title('欧拉法不同步长下输出响应曲线和真实响应曲线的比较');xlabel('t/time');ylabel('y/gain');legend('Euler 0.1步长','Euler 0.05步长','Euler 0.5步长','真实输出');欧拉法、二阶龙库、四阶龙库比较h=0.5%欧拉法求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(221)x=1;y=x;h=0.5;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'c')hold onx=1;y=x;h=0.5;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'m')title('欧拉法求出的输出响应曲线和真实的输出响应曲线之间的比较'); xlabel('t/time');ylabel('y/gain');legend('Euler','真实输出');%二阶龙格库塔法求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(222)x=1;y=x;h=0.5;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1);x=x+h*(k1+k2)/2;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'b')hold onx=1;y=x;h=0.5;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'g')title('二阶龙格库塔法求出的输出响应曲线和真实输出响应曲线之间的比较'); xlabel('t/time');ylabel('y/gain');legend('RK2','真实输出');%四阶龙格库塔法求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(223)x=1;y=x;h=0.5;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1/2);k3=-(x+h*k2/2);k4=-(x+h*k3);x=x+h*(k1+k2*2+k3*2+k4)/6;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'*g')hold onx=1;y=x;h=0.5;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'r')title('四阶龙格库塔法求出的的输出响应曲线和真实输出响应曲线之间的比较'); xlabel('t/time');ylabel('y/gain');legend('RK4','真实输出');%欧拉法、二阶龙格库塔法和四阶龙格库塔法分别求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(224)x=1;y=x;h=0.5;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'k')hold onx=1;y=x;h=0.5;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1);x=x+h*(k1+k2)/2;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'b')hold onx=1;y=x;h=0.5;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1/2);k3=-(x+h*k2/2);k4=-(x+h*k3);x=x+h*(k1+k2*2+k3*2+k4)/6;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'g')hold onx=1;y=x;h=0.5;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'k')title('Euler、RK2、RK4下输出响应曲线和真实输出响应曲线之间的比较'); xlabel('t/time');ylabel('y/gain');legend('Euler','RK2','RK4','真实输出');h=0.1%欧拉法求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(221)x=1;y=x;h=0.1;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'c')hold onx=1;y=x;h=0.1;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'m')title('欧拉法求出的输出响应曲线和真实的输出响应曲线之间的比较'); xlabel('t/time');ylabel('y/gain');legend('Euler','真实输出');%二阶龙格库塔法求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(222)x=1;y=x;h=0.1;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1);x=x+h*(k1+k2)/2;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'b')hold onx=1;y=x;h=0.1;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'g')title('二阶龙格库塔法求出的输出响应曲线和真实输出响应曲线之间的比较'); xlabel('t/time');ylabel('y/gain');legend('RK2','真实输出');%四阶龙格库塔法求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(223)x=1;y=x;h=0.1;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1/2);k3=-(x+h*k2/2);k4=-(x+h*k3);x=x+h*(k1+k2*2+k3*2+k4)/6;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'*g')hold onx=1;y=x;h=0.1;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'r')title('四阶龙格库塔法求出的的输出响应曲线和真实输出响应曲线之间的比较');xlabel('t/time');ylabel('y/gain');legend('RK4','真实输出');%欧拉法、二阶龙格库塔法和四阶龙格库塔法分别求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(224)x=1;y=x;h=0.1;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'k')hold onx=1;y=x;h=0.1;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1);x=x+h*(k1+k2)/2;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'b')hold onx=1;y=x;h=0.1;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1/2);k3=-(x+h*k2/2);k4=-(x+h*k3);x=x+h*(k1+k2*2+k3*2+k4)/6;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'g')hold onx=1;y=x;h=0.1;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'k')title('Euler、RK2、RK4下输出响应曲线和真实输出响应曲线之间的比较'); xlabel('t/time');ylabel('y/gain');legend('Euler','RK2','RK4','真实输出');h=0.05%欧拉法求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(221)x=1;y=x;h=0.05;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'c')hold onx=1;y=x;h=0.05;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'m')title('欧拉法求出的输出响应曲线和真实的输出响应曲线之间的比较'); xlabel('t/time');ylabel('y/gain');legend('Euler','真实输出');%二阶龙格库塔法求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(222)x=1;y=x;h=0.05;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1);x=x+h*(k1+k2)/2;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'b')hold onx=1;y=x;h=0.05;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'g')title('二阶龙格库塔法求出的输出响应曲线和真实输出响应曲线之间的比较'); xlabel('t/time');ylabel('y/gain');legend('RK2','真实输出');%四阶龙格库塔法求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(223)x=1;y=x;h=0.05;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1/2);k3=-(x+h*k2/2);k4=-(x+h*k3);x=x+h*(k1+k2*2+k3*2+k4)/6;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'*g')hold onx=1;y=x;h=0.05;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'r')title('四阶龙格库塔法求出的的输出响应曲线和真实输出响应曲线之间的比较');xlabel('t/time');ylabel('y/gain');legend('RK4','真实输出');%欧拉法、二阶龙格库塔法和四阶龙格库塔法分别求出的输出响应曲线和真实的输出响应曲线之间的比较subplot(224)x=1;y=x;h=0.05;N=round(1/h);%取整t=0;T=0;for i=1:Nk=-x;x=x+h*k;y=[y,x];T=[T,t+h];%画图的时候都保存了初始值t=t+h;endplot(T,y,'k')hold onx=1;y=x;h=0.05;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1);x=x+h*(k1+k2)/2;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'b')hold onx=1;y=x;h=0.05;N=round(1/h);%取整t=0;T=0;for i=1:Nk1=-x;k2=-(x+h*k1/2);k3=-(x+h*k2/2);k4=-(x+h*k3);x=x+h*(k1+k2*2+k3*2+k4)/6;y=[y,x];T=[T,t+h];t=t+h;endplot(T,y,'g')hold onx=1;y=x;h=0.05;N=round(1/h);t=h;T=0;for i=1:Nx=exp(-t);y=[y,x];T=[T,t];t=t+h;endplot(T,y,'k')title('Euler、RK2、RK4下输出响应曲线和真实输出响应曲线之间的比较'); xlabel('t/time');ylabel('y/gain');legend('Euler','RK2','RK4','真实输出');。

微分方程数值解4.1当常微分方程能解析求解时，可利用Matlab符号工具箱中的功能找到精确解. 见下例求解方程,,,. 键入： yyy，，,20syms x y %定义符号变量diff_equ= ‘D2y+2*Dy-y=0’; %D2y表示,,,Dy= y,yy=dsolve (diff_equ, ‘x’) %定义x为自变量 y=cl*exp ((2^(1/2)-1)*x+c2*exp (-(2^(1/2)+1)*x)%表达式中含c1与c2，表示通解.%初始条件为y (0)=0,,y(0)=1时，按如下方式调用 y=dsolve (diff_equ,‘y (0)=0’, ‘Dy (0)=1’, ‘x’) y=1/4*2^(1/2)*exp ((2^(1/2)-1)*x)—1/4*2^(1/2)*exp (-(2^(1/2)+1)*x)%画出函数y=y (x)的图形ezplot (y,[-2,2])图形具体形式请上机试之.在方程无法获得解析解的情况下，可方便地获得数值解. 下面的例子说明用Matlab求数值解的方法及应注意的问题.[例1] 求解范德堡（vander pol）方程2dxdx2,,，,,(1)0xx 2dtdt求解高阶方程，必须等价地变换为一阶微分方程组，对本例，通过定义两个新的变量，实现这一变换yxydxdt1,2/,, 则令 dydty1/2,2dydtyyy2/(11)*21,,,,编写求解程序分为两部分，第一部分为待求解的方程，存盘的文件名为,待求解方程的函数名.m,，第二部分为求解主程序，本例中取名为main1.m.首先编写待求解方程的文件. 文件存盘名为“vdpol.m”. M,function yprime=vdpol(,)tyyprime (1)=y (2);; (1(1)^2)*(2)(1),,yyymu=2 yprime=[yprime (1);yprime (2)]; yprime (2)=mu*说明函数yprime=vdpol中. 定义为自变量，的形式取决于求解方程的阶数，本(,)tyyt 例中，,为解向量，为导数向量. yprime， y(2)yyyy,[(1),(2)],(1)(1)(1),y yprime,,，函数返回vander pol方程的导数列向量. 因为所求结果为方程数值解，(2)(1),y所以各向量维数只有在主程序求解时定下精度后才能确定.主程序定名为main1.m，你可用你所喜欢的其它名子，但vdpol.m除外. clear functions%调试程序时，放置这一语句是必要的. 它清除前边已编译的存在于内存中的废弃程序[]=ode23 (‘vdpol’,[0,30],[1,0]); ty,y1=y (:,1); %解曲线.y2=y (:,2); %解曲线的导数.polt ( ‘_ _’) tyty,1,,2,说明龙格_库塔的2阶与4阶改进型求解公式的实现，其指令分别为：[]=ode23 (‘f’,tsx,0,options) tx,[]=ode45 (‘f’,tsx,0,options) tx,其中可由系统依据精度要求自动设定，亦可由使用者依据实际需要自己确定，分别说明之. ts（1）若令tstttf,[0,1,,]，则输出在指定时刻tttf0,1,,给出，当tstktf,0::时，输出在区间[0,]ttf的等分点上给出，为步长. k（2）若tsttft,[0,],0为自变量初值，tf为终值，此时，options决定自变量的维数，t中的时间点不是等间隔的，这是为了保证所需的相对精度，积分算法改变了步长. 用于t,3,6设定误差限的参数options可缺省，此时系统设定相对误差为，绝对误差为，若1010自行设定误差限，可用如下语句：options=odeset (‘reltol’,, ‘abstol’,) rtat这里的与分别为设定的相对与绝对误差. rtat须注意的是无论用哪种方法确定ttf0,的取值方式，必须由使用者确定且应与相匹配. x0t,y0[1,0],ts,[0,30]y(0)1,y(0)0,为初始条件，本例中，因为，这意味着解曲线，，x0一般说，当解nnn个未知函数的方程组时，为维向量，共含有个初始条件. x0两个输出参数是列向量xx与矩阵，它们具有相同的行数，而矩阵的列数等于方程t组的个数，本例中y(:,1)y(:,2)的列数为2，其中，为自变量上各点函数值，为上各ytt点导数值.最后，提请读者注意的是：ode45也不总是比ode23好，在很多时候，低阶算法更有效，有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍. 有些参考书提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法. 4.2 -设有一阶方程与初始条件,yfxy,(,), （4.1） ,yxy(),00,其中适当光滑，关于满足Lipschitz条件，即存在使 fxy(,)LyfxyfxyLyy(,)(,),,,1212则（4.1）式的解存在且惟一.关于yyx,()的解析解一般难以求到或根本无解析解，因而，实际问题中，通常，采用差分的方法. 在一系列离散点xxx,,,,yyy,,,,上寻求其数值近似解. 12n12n相邻两个节点间的间距xxnhn,，,,1,2,hxx,,称为步长，一般地取等步长，则hn0nnn，11、欧拉方法在区间[,]xx上用差商 nn，1yxyx()(),nn，1 h代替（4.1）式中,[,]xxxxxy，对fxy(,)中在上取值还是，而形成向前欧拉公式nn，1nn，1与向后欧拉公式.（1）向前欧拉公式xfxy(,)取左端点，得如下公式 nyxyxhfxyx()()(,()),，（4.2） nnnn，1从yxy(),x点出发，由初值代入（4.2）求得 000yyhfxy,，(,) （4.3） 1000反复利用（4.2），有yyhfxyn,，,(,) 0,1,2, （4.4） nnnn，1其几何意义如图4.1所示.y 图中yyx,()为方程（4.1）的精确 P P43P 2解曲线，其上任意点(,)xy处切线斜率为误差 P 1yyx,() 32Pxy(,)fxy(,). 从初值点出发，用该 P000 0y 0点斜率fxy(,)xx,作一直线段，在 001yyx,() yx() 3处得到Pxy(,)y，由（4.2）式确定， 1111y 3再从Pxy(,)fxy(,)出发，以为斜率 11111作直线段，在xx,Pxy(,)处得到， 2222xxxxx x O03124PPP， 012作为积分曲线yx()的近似，用图4.1 yyx,()n，1这一过程继续下去，形成折线表示在xy处的精确值，为解的近似值，不难得到 n，1n，12h32,,()()()()yxyyxOhOh,,，, nnn，，112P，1这一误差称为局部截断误差. 若一种算法局部截断误差为Oh()，则称该算法具有阶P精度，所以向前欧拉公式具有1阶精度.（2）向后欧拉公式若[,]xxxx中取中的，则有如下公式： fxy(,)nn，1n，1yyhfxyn,，,(,) 0,1,2, （4.5） nnnn，，，111称式（4.5）为向后欧拉公式，因为此式中y未知，故称其为隐式公式，无法用其直n，1接计算y，一般用向前欧拉公式产生初值. n，1(0)yyhfxyn,，,(,) 0,1,2, 11nnnn，，再按下式迭代(1)()kk，yyhfxykn,，,,(,),0,1,,0,1, nnnn，，，111其误差估计如下2h32,,()()()()yxyyxohoh,,，, nnn，，112精度亦为1阶，将向前欧拉公式（4.4）与向后欧拉公式（4.5）及它们的误差的几何说明作一对比，是十分有益的，见图4.2.y 为讨论局部截断误差，在图4.2中设点APxy(,)落在积分曲线yyx,()上，按式 nnnyyx,() （4.4）及式（4.5）分别得 ,P点为与， ABn，1 B且P AB,yyx,()点一定在积分曲线上相应 n点的上、下两边,所以将式（4.4）与（4.5） , 平均之，一定能得到更好的结果.xxx （3）梯形公式 nn，1 将向前与向后欧拉公式加以平均得到所图4.2 谓梯形公式hyyfxyfxyn,，，,[(,)(,)] 0,1,2, （4.6） nnnnnn，，，11123其局部截断误差为Oh()，具有2阶精度.（4）改进的欧拉公式为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步yyhfxy,，(,)nnnn，1h （4.7） yyfxyfxyn,，，,[(,)(,)], 0,1,2,nnnnn，，11n，12或写为h,yykk,，，()nn，112,2,1nn （4.8） n,0,1,2,kfxy,(,),,211nn，kfxyhk,，(,),,最后指出，上述欧拉方法可推广至微分方程组，如,yfxyz,(,,),,,zgxyz,(,,), ,yxy(),00,,zxz(),,00向前欧拉公式为yyhfxyz,，(,,),nnnnn，1 n,0,1,2, ,zzhgxyz,，(,,),nnnnn，12、龙格_库塔方法由微分中值定理,[()()]/(),01yxyxhyxh,,，,,,, nnn，1又因为,,yxhfxhyxh()(,())，,，，,,,yfxy,(,)，所以 nnn从而有yxyxhfxhyxh()()(),(),，，，,, （4.9） nnnn，1令[,]xx，称其为区间上的平均斜率，由（4.9）可知，给kfxhyxh,，，(,()),,nn，1nn出一种平均斜率，可相应导出一种算法. 向前欧拉公式中，精度低. 改进欧kfxy,(,)nn1拉公式中取[,]xxkfxyfxy,，((,)(,))，精度提高，下面，我们在区间内nn，1nnn，1n，12多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体给出，只开列具体结果.（1）2阶龙格_库塔公式yyhkk,，，(),,,nn，11122,kfxy,(,) （4.10） 1nn,,21nnkfxahyhka,，，,,(,),0,1,,,1,其中,,,，,,,1,,1a，由于4个未知数只有3个方程，所以解不惟一，若令1222a1，即得改进的欧拉公式，具有2阶精度. ,,,，,,,,1a122（3）4阶龙格_库塔公式只给出精典格式中最常用的一种.h,yykkkk,，，，，(22)nn，11234,6,kfxy,(,),1nn,hh, （4.11） kfxyk,，，(,),21nn22,hh,kfxyk,，，(,)32nn,22,kfxhyhk,，，(,),43nn,其计算精度为4阶4.31、模型与问题[例2] 单摆运动图4.3中一根长的细线，一端固定，另一端悬挂质量为 lm的小球，在重力作用下，小球处于竖直的平衡位置. 现使小球偏离平衡位置一个小的角度，然后使其自由运动，在不 ,考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动.为平衡位置，在小球摆动过程中，当与平衡位置夹 ,,0角为,mgsin,时，小球所受重力在其动运轨迹的分量为 , , l（负号表示力的方向使减少），由牛顿第二定律可得微分方 ,程,,mltmg,,()sin,, （4.12）设小球初始偏离角度为,，且初速为0，式（4.12）的初 0,始条件为,,,,(0),(0)0,, （4.13） 0 mg 当,不大时，，式（4.12）化为线性常系数微sin,,,0分方程图4.3g,,,,，,0 （4.14） l解得g,,()costt, （4.15） 0ll简谐运动的周期为. T,2,g现在的问题是：当,较大时，仍用近似，误差太大，式（4.12）又无解析解，,sin,0试用数值方法在,,30,10,,两种情况下求解，画出的图形，与近似解（4.15）,()t00比较，这里设. l,25cm[例3] 捕食与被捕食当鲨鱼捕食小鱼，简记为乙捕食甲，在时刻，小鱼的数量为，鲨鱼的数量为，xt()yt()t当甲独立生存时它的（相对）增长率与种群数量成正比，即有,r，为增长率，xtrxt()(),而乙的存在使甲的增长率r减少，设减少率与乙的数量成正比，而得微分方程,xtxtraytrxaxy()()(()),,,, （4.16）比例系数a反映捕食者掠取食饵的能力.乙离开甲无法生存，设乙独自存在时死亡率为,，ytdyt()(),,，甲为乙提供食物，d使乙的死亡率降低，而促其数量增长，这一作用与甲的数量成正比，于是yt()满足 d,ytydbxtdybxy()(()),,，,,，（4.17）比例系数反映甲对乙的供养能力，设若甲，乙的初始数量分别为 bxxyy(0);(0),, （4.18） 00则微分方程（4.16），（4.17）及初始条件（4.18）确定了甲，乙数量xt()、yt()随时间变化而演变的过程，但该方程无解析解，试用数值解讨论以下问题：（1）设rdabxy,,,,,,1,0.5,0.1,0.02,25,2，求方程（4.16），（4.17）在00条件（4.18）下的数值解，画出xtyt(),()的图形及相图(,)xy，观察解xtyt(),()的周期变化，近似确定解的周期和的最大、小值，近似计算在一个周期内的平均值. xy,xy,（2）从式（4.16）和（4.17）消去得到 dtdxxray(),, （4.19） dyydbx(),，解方程（4.19），得到的解即为相轨线，说明这是封闭曲线，即解确为周期函数.（3）将方程（4.17）改写为,1yxtd()[],，（4.20） by在一个周期内积分，得到xt()yt()一周期内的平均值，类似可得一周期内的平均值，将近似计算的结果与理论值比较.2、实验（1）方程（单摆问题）,,mlmg,,,,sin, ,,,,,(0),(0)0,,0,无解析解，为求其数值解，先将其化为等价的一阶方程组. 令,xx,,,,,，方程化为 12,,xx,12,g,,21,,xxsin ,l,102,,xx(0),(0)0,,,其中，gl,,9.8,25,,为（弧度）与（弧度）两种情况，具100.1745,300.5236,0体编程如下：先以danbai.m为文件名存放待解方程. 键入： function xdot =danbai (t,x) %x=[x (1),x (2)],g=9.8;1=25; %x (1)为解向量, x (2)是导数 xdot (1)=x (2); %xdot(1)=,,(1)= x,xdot (2)=-g/1*sin (x (1)); %xdot (2)=,, ,xdot=[xdot (1);xdot (2)]; %必须将导数向量以列向量形式给出.再以主程序(求数值解)调用待求方程，主程序用main2.m为文件名存盘，其代码如下：clear functions[t,x]=ode23 (‘danbai’,[0,10],[0.1745,0])%只计算,,100.1745,（弧度）的情形.01g,,()cos,twtw,,%对近似解，周期 T2,,01gw=sqrt (9.8/25);y=0.1745*cos (w*t);[t,x (:,1),y] %显示数据，无分号. hold on %欲在同一幅图中画近似解. plot (t,x (:,1), ‘b’) %画数值解，绿色plot (t,y, ‘g*’) %用*号，红色画近似解. Hold off（2）食饵_捕食者方程,xtrxaxy(),,,ytdybxy(),,，可化为如下形式，，xrayx,0，，，，,,, ,,,,,,0,，dbxy，，，，y，，初始条件xxyy(0),(0),,表示为 00xx(0)，，，，0 ,,,,,yy(0)，，，，0以shier.m存盘如下代码function xdot=shier (t,x)r=1;d=0.5;a=0.1;b=0.02;xdot=diag ([r-a*x (2),-d+b*x (1)])*x;，，xxx(1)，，，，%x=,,,xdot= ,,,,,,,xy(2)，，，，y，，以main3.m存盘如下代码.clear functionsts=0:0.1:15;x0=[25,2];[t,x]=ode45 (‘shier’,ts,x0);[t,x] %显示数据t,x,y plot (t,x)gird %加网格线gtext (‘x(t)’),gtext (‘y(t)’), %用点鼠标方式pause,figure (2) %将x(t),x(t)放至指定点 12plot (x (:,1),x (:,2)) %以x (1)与x (2)为坐标点画相图xlabel (‘x’),ylabel (‘y’)可以猜测，xt()xt()(,)xx与是周期函数，相图是封闭曲线，从数值解可近似定出2121周期为10.7,x2.0,x的最大和最小值分别为与的最大和最小值分别为和，99.328.42.012为求xx与在一个周期的平均值，只需键入： 12y1=x (1:108,1); %x周期为10.7. 1x1p=trapz (y1)*0.1/10.7, %trapz (y1)返回按 y2=x (1:108,2); %梯形法对y1的积107yiyi1()1(1)，，x2p=trapz (y2)*0.1/10.7, %分值,x ,i12i,可得xx,,25,10 12对方程dxxray(),,，,dbxray,dxdy,化为 dyydbx(),，xy积分得解dbxray,,()()xeyec, （4.21）即为原方程组的相轨线，其中c由初始条件确定. 为说明上述相轨线是封闭的，令dbxray,,fxxegyye();(),,设其最大值分别为fg,xy,，若满足 mm00fxffyg(),(),, 00mmdr则有,,xy,fgc,,（令，可解出），又当时，相轨线（4.21）fg,,0,0xy,,,00mm00ba有意义. 当fgc,0,,cfgPxy(,)时，相轨线退化为一个点，对时，相轨线如mmmm00图4.4（c），而图（a），（b）分别为fx()与gy()的图像，下面讨论如何由（a），（b）画（c）.fxf(),nm y gyg(),gv()nm fx() y2 gm fm Pxy(,)qq 00 1 yP q02 q3 y1x y x x xxxy xyxy012001221 （a）（b）（c）图4.4设cpgpf,,,(0)yy,，xx,fxp(),，若令，则有，由（a）知，，使mm012qxy(,)xxx,,qxy(,)fxfxp()(),,，且于是相轨线应过与，对11010222012 fxgxpggyg()()(),,,xxx,(,)fxp(),gyq(),，，由，令，又由（b）知，mm12 存在yyy,,qxy(,)yy,gygyq()(),,qxy(,)使，且，于是相轨线又过与10231121242yyyyy,(),,xqq,两点，所以对间每个，轨线总要通过纵坐标为的两点，于是1210212相轨线是一条封闭曲线.相轨线封闭等价于xtyt(),()是周期函数，记周期为，为求其在一个周期内的平均T值(,)xy，由1yxtd()(),， by两边在一个周期内积分有((0)())yyT,：T11ln()ln(0)yTydTd,，，xxtdt,,，,() ,,,0TTbbb，，同理ry, a从而xxyy,,,00即xy的平均值恰为相轨线中心点的坐标，提请读者注意的是：这里的与与xtyt(),()p00初始条件中的xy,不是一件事. 00注意到在生态学上的意义，上述结果表明，捕食者的数量与食饵的增长率成rdab,,,正比，与它掠取食饵的能力成反比，食饵的数量与捕食者死亡率成正比，与它供养捕食者的能力成反比.3、练习内容（1）编写改进欧拉公式求微分方程数值解的程序，并用其与ode23求下列微分方程数值解，对二者作出比较.a）22,yxyy,,,,(0)0或y(0)1,.222 b),,,xyxyxny，，,,()02,,,21 yxsin,n,,（Bessel方程，这里令，其精确解为）. yy()2,(),,,x222,c),,,yyxyy，,,,cos0,(0)1,(0)0.（2）倒圆锥形容器，上底面直径为1.2m，容器的高亦为1.2m，在锥尖的地方开有一直径为3cm的小孔，容器装满水后，下方小孔开启，由水利学知识可知当水面高度为时，h水从小孔中流出的速度为为重力加速度，若孔口收缩系数为0.6（即若一个面vghg,2,积单位的小孔向外出水时，水柱截面积为0.6），问水从小孔中流完需多少时间？2分钟时，水面高度是多少？（3）一只小船渡过宽为的河流，目标是起点正对着的另一岸上点，已知河水ABd流速vv与船在静水中的速度之比为. k12（a）建立小船航线的方程，求其解析解.（b）设dvv,,,100m,1m/s,2m/s，用数值解法求渡河所需时间，任意时刻小12 船的位置及航行曲线，作图并与解析解比较.（c）若流速v为0,0.5,2 (m/s)，结果将如何？ 1（4）研究种群竞争模型. 当甲、乙两个种群各自生存时，数量演变服从下面规律xyxtrxytry()(1),()(1),,,, 12nn12其中，rr,nn,xtyt(),()分别为时刻甲，乙两个种群的数量，为其固有增长率，为它t1212们的最大容量，而当这两个种群在同一环境中生存时，由于乙消耗有限资源而对甲的增长产生影响，将甲的方程修改为xyxrxs,,,(1) （4.22） 11nn12这里s的含意是：对于供养甲的资源而言，单位数量乙（相对n）的消耗率为单位数12量甲（相对n）消耗的s倍，类似地，甲的存在亦影响乙的增长，乙的方程应改为 11xyytrys()(1),,, （4.23） 22nn12给定种群的初始值为xxyy(0),(0),, （4.24） 00及参数rrssnn,,,,,后，方程（4.22）与（4.23）确定了两种群的变化规律，因其解析121212解不存在，试用数值解法研究以下问题：（a）设rrnnssxy,,,,,,,,1,100,0.5,2,10，计算，画出xtyt(),()12121200 它们的图形及相图(,)xy，说明时间充分大以后xtyt(),()的变化趋势（人们今天看到的t已经是自然界长期演变的结果）.（b）改变rrxy,,,，但s与s不变（保持ss,,1,1），分析所得结果，若12001212ss,,,,1.5(1),0.7(1)，再分析结果，由此你得到什么结论，请用各参数生态学上的12含义作出解释.（c）试验当ss,,,,0.8(1),0.7(1)ss,,,,1.5(1),1.7(1)时会有什么结果；当1212时又会出现什么结果，能解释这些结果吗？。

matlab求解二阶微分方程和一阶微分方程的平方代码一、介绍在数学和工程领域，微分方程是一种常见的数学工具，用于描述一些变化的规律和规律性。

其中，二阶微分方程和一阶微分方程是极为重要的类型，广泛应用于控制系统、信号处理、机器学习等领域。

而在实际的工程问题中，利用 MATLAB 求解二阶微分方程和一阶微分方程是一种常见且高效的方法。

本文将重点探讨如何使用 MATLAB 求解这两类微分方程，并给出相应的代码示例。

二、二阶微分方程的求解对于形如 y'' = f(x, y, y') 的二阶微分方程，可以通过 MATLAB 的ode45 函数进行求解。

ode45 是 MATLAB 中最常用的数值求解微分方程的函数，其基本用法如下：```matlabfunction dy = myODE(x, y)dy = zeros(2, 1);dy(1) = y(2);dy(2) = -0.1 * y(2) - sin(y(1));end[t, y] = ode45(@myODE, [0, 20], [0, 1]);plot(t, y(:, 1));```在该示例中，我们定义了一个名为 myODE 的函数来描述 y'' = -0.1y' - sin(y)，然后使用 ode45 函数对其进行数值求解，并绘制出了解 y 关于 x 的图像。

三、一阶微分方程的求解对于形如 dy/dx = f(x, y) 的一阶微分方程，同样可以利用 ode45 函数进行求解。

其基本用法如下：```matlabfunction dy = myODE(x, y)dy = -2*x*y;end[t, y] = ode45(@myODE, [0, 1], 1);plot(t, y);```在该示例中，我们定义了一个名为 myODE 的函数来描述 dy/dx = -2xy，然后同样使用 ode45 函数对其进行数值求解，并绘制出了解 y 关于 x 的图像。

四阶龙格库塔方法求解n自由度二阶微分方程在数值计算中，龙格-库塔方法（Runge-Kutta method）是一种常用的求解常微分方程（ODE）的数值方法。

它是一种多步法（multiple-step method），通过使用不同的近似值来迭代计算ODE 的解。

其中，四阶龙格-库塔方法是最常用的龙格-库塔方法之一，也是一种高阶方法，可用于求解n自由度二阶微分方程。

下面将介绍四阶龙格-库塔方法的原理、步骤和应用。

首先，我们来回顾一下二阶微分方程的一般形式：y''(y)=y(y,y(y),y'(y))这里，y(y,y(y),y'(y))是已知的函数，y(y)是我们要求解的未知函数。

为了求解这个二阶微分方程，我们需要将其转化为一个一阶微分方程的求解问题。

令y(y)=y(y)y(y)=y'(y)这样，我们可以将原始的二阶微分方程转化为如下的一阶微分方程组：y'(y)=y(y)y'(y)=y(y,y(y),y(y))现在，我们可以利用四阶龙格-库塔方法来求解这个一阶微分方程组。

四阶龙格-库塔方法基于泰勒展开的思想，通过使用一系列的近似值来逼近方程的解。

该方法的一般形式为：y1=ℎy(yy,yy)y2=ℎy(yy+ℎ/2,yy+y1/2)y3=ℎy(yy+ℎ/2,yy+y2/2)y4=ℎy(yy+ℎ,yy+y3)yy+1=yy+1/6(y1+2y2+2y3+y4)yy+1=yy+ℎ其中，yy是当前时间步，yy+1是下一个时间步，yy是当前位置的近似解，yy+1是下一个位置的近似解，ℎ是时间步长，y1、y2、y3和y4是中间的近似值。

四阶龙格-库塔方法的步骤如下：Step 1:给定初始条件我们需要提供初始条件，包括初始位置y0和初始速度y0，以及时间步长ℎ、求解的时间区间[y0,yy]和离散时间步数y。

y0=y(y0)y0=y'(y0)Step 2:进行迭代计算从n=0开始，进行y步的迭代计算，其中y=(yy−y0)/ℎ。

1. matlab 新建.m文件，编写龙格-库塔法求解函数function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）n=floor((b-a)/h); %求步数x(1)=a;%时间起点y(:,1)=y0;%赋初值，可以是向量，但是要注意维数for ii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii));k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);k4=ufunc(x(ii)+h,y(:,ii)+h*k3);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龙格库塔方法进行数值求解end2.另外再新建一个.,m文件，定义要求解的常微分方程函数function dx=fun1(t,x)dx =zeros(2,1);%初始化列向量dx(1) =0.08574*x(2)-1.8874*x(1)-20.17;dx(2) =1.8874*x(1)-0.08574*x(2)+115.16;3,再新建一个.m文件，利用龙格-库塔方法求解常微分方程[T1,F1]=runge_kutta1(@fun1,[46.30 1296 ],1,0,20); %求解步骤2定义的fun1常微分方程，@fun1是调用其函数指针，从0到20，间隔为1subplot(122)plot(T1,F1)%自编的龙格库塔函数效果title('自编的龙格库塔函数')grid运行步骤3文件即可得到结果，F1为估测值或者可以调用matlab自带函数ode45求解方法如下：[T,F] = ode45(@fun1,0:1:20,[17.64 37800 ]); subplot(121)plot(T,F)%Matlab自带的ode45函数效果title('ode45函数效果')grid。