Matlab微分方程数值解（专题七）_使用matlab自带的龙格-库塔方法求解下列函数y(x0)及其导数在区间[o 10]上的函

前言

对于某些非线性微分方程是没有函数解的，只能用数值解法求解，因为本文章介绍两种常规方法求数值解。

Runge-Kutta

[t,y] = ode23(odefun,tspan,y0)

此 MATLAB 函数（其中 tspan = [t0 tf]）求微分方程组 y'=f(t,y) 从 t0 到 tf 的积分，初始条件为 y0。解数组 y 中的每一行都与列向量 t 中返回的值相对应。

顾名思义就是该命令针对的是一阶y'=f(t,y)函数的形式。

积分起始值是t0,终止值是tf

y0是初始条件（或者可以说是初始状态列向量）

odefun为要求解的函数，需注意是y'=f(tf)的形式。编写形式可以用M函数，也可以用函数句炳形式。

例如求y'=5+（3y-4x)^2的(0≤x≤10)数值解，满足初始条件y(0)=0.5;

函数句柄形式

clear
ff=@(x,y)5+(3*y-4*x)^2;
[x,y]=ode23(ff,[0,10],0.5)

函数句柄不需要单引号，m函数文件需要 

M文件形式

建立c.m文件，记得保存

function ff=c(x,y)%yy为y'
ff=5+(3*y-4*x)^2;

命令行输入

clear
[x,y]=ode23('c',[0,10],0.5)

clear
f1=@(t,y)(y^2-t-2)/(4*(t+1));
[t,y]=ode23(f1,[0 10],2);
y1=sqrt(t+1)+1;%t是x轴上的坐标
A=[y,y1]%输出龙格数值解和解析解进行比较。
plot(t,y,'-r')
hold on 
plot(t,y1,'b--')
legend('龙格数值解','解析解')

 Van Der Pol

求著名的vanderpol方程的数值解，并绘制其时间对应的曲线和状态轨迹图。

由于ode法只能用作一阶，所以我们要把他看出一阶来做，并变成两个方程组形式

故设y1=x',y2=x;

所以y2'=y1

所以y1'=-x-x'(x^2-1)=y1(1-y2^2)-y2

所以我们编写M文件如下

function F=c(t,y)%yy为y'
F=ones(2,1);%要把F设为列向量。F=zeros(2,1)也可以
F(1)=(1-y(2)^2)*y(1)-y(2);
F(2)=y(1);

命令行输入

clear
[t,y]=ode23('c',[0 20],[0,0.25]);%t的范围0到20，两个y的初始值分别为0和0.25
plot(t,y(:,1),':b',t,y(:,2),'-r')%时间相应曲线
legend('速度（y1）','位移（y2）')
figure(2);
plot(y(:,1),y(:,2))%状态轨迹图

 Euler

向前欧拉法又称显式欧拉公式

编写欧拉公式函数文件如下

function P=Euler(x0,y0,b,h)
%y0为初始值
n=(b-x0)/h;%x0为初值，b为最右边的端点（x的最大值）
X=zeros(round(n),1);%取整是为了防止n作为小数出现。取整命令round
Y=zeros(round(n),1);
k=1;
X(k)=x0;
Y(k)=y0;
for k=1:n
    X(k+1)=X(k)+h;
    Y(k+1)=Y(k)+h*(X(k)-Y(k));
    k=k+1;
end
 
%前面是编写好的程序，下面求y=e^x-x
%满足初始条件y(0)=1的数值解。步长取0.1，x的最大值取5
y=exp(X)-1
plot(X,y(:,1),'mp')%数值解
%也可以plot(x,y,'mp'),上面是为了防止y出现矩阵形式，所以要确保拿到第一列准确数值。
legend('欧拉数值解')

接下来在命令行窗口输入以下，并与原函数 精确解进行图形比较`

oula(0,1,5,0.1)%因为我的文件名是oula，所以是这样子的格式，大家要根据自己的文件名来调用程序
hold on
x=0:0.1:5;
a=exp(x)-x%去掉分号是为了和上面输出的y值作为比较。
plot(x,a,'--g')

大家可以尝试把 步长0.1改成0.01，得到的结果可能会更精确。