Matlab 解常微分方程的初值问题
 
题目：Matlab 解常微分方程的初值问题
 
设计目的：
 
1、熟练掌握Matlab的基本编程方法，及其编程风格。
 
2、熟练掌握Matlab常用函数的使用。
 
3、与本专业相关知识相结合，掌握其在程序开发中的应用方法
 
以及和word、C语言等接口方法。
 
4、通过计算机数值求解的方式来加深微分方程解的理解。
 
5、熟悉初等方法可获得解析解之外的数值近似解的求解方法，提
 
高对差分格式的认识和离散化分析问题的技巧，加深对理论课程的学习和理解，为数学专业和信息与计算科学专业其他后继课程的学习打好基础。
 
设计内容：
 
已知一个三阶微分方程:，利用matlab软件求这个三阶微分方程在初值 下的解。
 
原三阶微分方程可化为：
 
令 则原三阶微分方程可化为微分方程组 在初值 下的解。
 
程序流程：
 
程序代码：
 
%编写函数文件rigid.m
 
function dy = rigid(t,y)
 
dy = zeros(3,1); % a column vector
 
dy(1) = y(2) ;
 
dy(2) =y(3);
 
dy(3) = 2*(1-y(1)^2)*y(3)-y(1)*y(2);
 
%调用函数ode45求解，时间区间为[0,10]
 
[t,Y] = ode45(@rigid,[0 10],[1 0 -1])
 
t =
 
0
 
0.0001
 
0.0001
 
0.0002
 
0.0002
 
0.0005
 
0.0007
 
………………
 
0.9383
 
1.0665
 
1.1947
 
1.2918
 
1.3889
 
1.4860
 
………………
 
6.2916
 
6.2922
 
6.2928
 
6.2934
 
6.2940
 
6.2947
 
6.
 
由于数据太多，这里只列举部分
 
%绘制解的曲线
 
plot(t,Y(:,1),'-',t,Y(:,2),'-.',t,Y(:,3),'.')
 
%给图形加标注
 
title('Solution of Rigid Equation')
 
xlabel('time T')
 
ylabel('solution Y')
 
legend('Y1','Y2','Y3')
 
设计结果：
 
未加图形标注时的图
 
加了图形标注后的图
 
结果分析：
 
输出结果[T,Y]中T为时间点组成的向量。Y为对应于T中时间点的y(1)、y(2)和y(3)的值。此次利用matlab数值方法来求解微分方程主要是把求解的时间划分成有限步，对应于每一步将计算出一个解，如果求得的解不满足误差限制，则减少步长，再求解。如此重复，直到满足误差限为止。
 
课设总结：
 
最初拿到题目后，开始读题，知道题目所表达的意思及我们所要完成的目的，达到的效果后，便开始了做题。
 
首先是找一个三阶微分方程。然后把它化简为标准形式，再利用matlab软件求解。其中，在利用matlab求解时遇到一些问题，比如画图时调用已经编号的rigid函数时的调用格式不正确，还有就是给图形家标注时程序的引号没有切换成英文输入法状态下的等一些问题。但是经过多次调试，检查，修改后程序运行总算成功。
 
通过此次课程设计，我已经基本掌握Matlab的基本编程方法，及其编程风格。且能较熟练掌握Matlab常用函数的使用。此次设计后我的收获不少。
 
参考文献：
 
[1] 张圣勤编 MATLAB7.0 机械工业出版社
 
[2]周义仓 靳祯 秦军林编 常微分方程极其应用 科学出版社
 
[3]韩明 王家宝 李林编 数学实验(matlab版) 同济大学出版社
 
[4]汪晓银 皱庭荣编 数学软件与数学实验 科学出版社
 
2、把这个三阶微分方程化为形如 的标准形式
 
1、已知一个三阶微分方程
 
3、编写函数文件rigid.m
 
4、调用函数文件rigid.m，利用ode45求解
 
5、绘制解的曲线

欧拉方法yn+1 = yn +hf(xn,yn) 是一种显式算法，计算量小，但精度很低。梯形方法yn+1 = yn + h / 2 *[f(xn,yn) + f(xn+1,yn+1)] 虽然提高了精度，但它是一种隐式算法，需要借助于迭代过程求解，计算量大。
 
为此，我们不能想到，综合应用这两种方法，先用欧拉格式求得一个初步的近似值yn+1′，称之为预报值；预报值y0的精度不高，我们用它替代梯形格式右端的yn+1再直接计算，的校正值yn+1.这样建立的预报-校正系统称作改进的欧拉格式：
 
 其可表示为如下嵌套格式：
 

 也可表示为下列平均化格式（便于编程实现）:
 

 程序框图如下：
 
 
例：用改进欧拉方法求解以下常微分方程的初值问题。
 
 其中取步长h=0.1
 
运行示例：
 
 
源码：
 
#include<iostream>
#include<iomanip>
double f(double x, double y);   
double f1(double x);

using namespace std;

int main(void)
{
	double x0, y0, h;   //x0,y0为初值;h为步长
	double x1, y1, yp, yc;
	int N;   //N为步数

	cout << "请输入初值:";
	cin >> x0 >> y0;

	cout << "请输入步长:";
	cin >> h;

	cout << "请输入步数:";
	cin >> N;

	//输出提示信息
	int i = 1;
	cout << "\t" << setw(10) << "xn" << "\t" << setw(10) << "yn" << "\t" << setw(10) << "y(xn)" << endl;

	for (i; i <= N; i++)
	{
		x1 = x0 + h;   //离散点
		yp = y0 + h * f(x0, y0);   //根据欧拉格式求出的值，预报
		//校正
		yc = y0 + h * f(x1, yp);    
		y1 = (yp + yc) / 2;
		//输出本次步进后的离散点数据
		cout << i << "\t" << setw(10) << x1 << "\t" << setw(10) << y1 << "\t" << setw(10) << f1(x1) << endl;
		//把新的值赋给旧值，开启下一趟循环
		x0 = x1;
		y0 = y1;
	}

	return 0;
}

double f(double x, double y)   //欧拉格式中的x，y的函数关系式，即f(xn,yn)
{
	double result;

	result = y - 2 * x / y;

	return result;
}

double f1(double x)   //实际函数解析式
{
	double result;

	result = sqrt(1 + 2 * x);

	return result;
}

 
求解初值问题
 
简介
前期准备
欧拉法
改进的欧拉法
龙格-库塔法
标准四阶显式Kutta公式
三级三阶显式公式
四级四阶显式Kutta公式
四级四阶显式Gill公式
  
示例
MATLAB代码
结果
 

 
简介
 
通过求解简单的初值问题：
 $$\{\begin{array}{llccccc}\frac{du}{dx}=f(x,u)& & & & & & (1)\\ u(x_0)=u_0& & & & & & (2)\end{array}$$
 引入欧拉法、改进的欧拉法、龙格-库塔法等。
 
前期准备
 
数值解法的基本思想就是先对$x$和$u(x)$在区间$[x_0,\infty )$上进行离散化，然后构造递推公式，再进一步得到$u(x)$在这些位置的近似取值。
 
取定步长$h$，令$x_n=x_0+nh(n=\pm 1,\pm 2,\cdots )$
得到离散的位置：$x_1,x_2,\cdots ,x_n,\cdots$
$u(x)$在这些点精确取值为：$u(x_1),u(x_2),\cdots ,u(x_n),\cdots$
利用数值解法得到的这些点的近似取值，记为$u_1,u_2,\cdots ,u_n,\cdots$
 
欧拉法
 
欧拉法的核心就是将导数近似为差商。
 将导数近似为向前差商，则有：
 $${\frac{du}{dx}|}_{x=x_n}\approx \frac{u(x_{n+1})-u(x_n)}{h}$$
 代入(1)式，有：
 $$u(x_{n+1})=y(x_n)+hf(x_n,u(x_n))$$
 用$u_{n+1}$和$u_n$代替$u(x_{n+1})$和$u(x_n)$，得：
 $$u_{n+1}=u_n+hf(x_n,u_n)$$
 因此，若知道$u_0$我们就可以递归出$u_1,u_2,\cdots$
 如果将导数近似为向后差商:
 $${\frac{du}{dx}|}_{x=x_n}\approx \frac{u(x_n)-u(x_{n-1})}{h}$$
 类似的，就可以得到：
 $$u_{n-1}=u_n-hf(x_n,u_n)$$
 这样，若知道$u_0$我们就可以递归出$u_{-1},u_{-2},\cdots$
 
改进的欧拉法
 
对$(1)$式在$[x_n,x_{n+1}]$上积分，可得：
 $$u(x_{n+1})=u(x_n)+{\int }_{x_n}^{x_{n+1}}f(x,u)dx$$
 其中，$n=0,1,\cdots$。用不同方式来近似上式的积分运算，就会得到不同的递推公式。若使用左端点计算矩形面积并取近似：
 $${\int }_{x_n}^{x_{n+1}}f(x,u)dx\approx hf(x_{n+1},u(x_{n+1}))$$
 代入上式得：
 $$u_{n+1}=u_n+hf(x_n,u_n)$$
 若使用梯形的面积做近似：
 $${\int }_{x_n}^{x_{n+1}}f(x,y)dx\approx \frac{h}{2}[f(x_n,u(x_n))+f(x_{n+1},u(x_{n+1}))]$$
 得到：
 $$u_{n+1}=u_n+\frac{h}{2}[f(x_n,u_n)+f(x_{n+1},u_{n+1})]$$
 欧拉法虽然精度偏低，但它是显式的，可直接得到结果。而梯形公式是隐式的，虽然精度较高，却无法通过一步计算得到结果，若用迭代法计算，运算量较大。综合这两种方法，可以相得益彰：先用显式格式却低精度的欧拉法计算得到一个粗略的预测值${\bar{u}}_{n+1}$，再将这个预测值代入梯形公式进行修正，得到较高精度的结果$u_{n+1}$。
 $$\{\begin{array}{l}{\bar{u}}_{n+1}=u_n+hf(x_n,u_n)\\ u_{n+1}=u_n+\frac{h}{2}[f(x_n,u_n)+f(x_{n+1},{\bar{u}}_{n+1})]\end{array}$$
 
龙格-库塔法
 
将以上两种方法分别写成如下形式：
 $$\{\begin{array}{l}{u}_{n+1}=u_n+h{K}_1\\ K_1=f(x_n,u_n)\end{array}$$
 $$\{\begin{array}{l}{u}_{n+1}=u_n+\frac{h}{2}(K_1+K_2)\\ K_1=f(x_n,u_n)\\ K_2=f(x_n+h,u_n+K_1)\end{array}$$
 上述方法都是通过$f(x,u)$在不同位置的线性组合来计算$u_{n+1}$的值，所考虑的位置越多，精度也越高。类似的，就得到龙格-库塔法的思想：如果用$f(x,u)$在更多位置的线性组合来构造递推公式，将会得到更高的精度。
 这样，递推公式将有如下形式：
 $$\{\begin{array}{l}{u}_{n+1}=u_n+h\sum _{i=1}^r{R}_i{K}_i\\ K_1=f(x_n,u_n)\\ K_i=f(x_n+a_{i}h,u_n+\sum _{j=1}^{i-1}{b}_{ij}{K}_j),i=2,3,\cdots ,r\end{array}$$
 其中，$R_i,a_i,b_{ij}$为待定常数。(利用$Taylor$展开就可以确定待定系数)
 
标准四阶显式Kutta公式
 
$$\{\begin{array}{l}{y}_{n+1}=y_n+\frac{h}{6}(K_1+4K_2+K_3),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{K}_1),\\ K_3=f(x_n+h,y_n-h{K}_1+2h{K}_2);\end{array}$$
 
三级三阶显式公式
 
$$\{\begin{array}{l}{y}_{n+1}=y_n+\frac{h}{4}(K_1+3K_3),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{3}h,y_n+\frac{1}{3}h{K}_1),\\ K_3=f(x_n+\frac{2}{3}h,y_n+\frac{2}{3}h{K}_2);\end{array}$$
 
四级四阶显式Kutta公式
 
$$\{\begin{array}{l}{y}_{n+1}=y_n+\frac{h}{8}(K_1+3K_2+3K_3+K_4),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{3}h,y_n+\frac{1}{3}h{K}_1),\\ K_3=f(x_n+\frac{2}{3}h,y_n-\frac{2}{3}h{K}_1+h{K}_2),\\ K_4=f(x_n+h,y_n+h{K}_1-h{K}_2+h{K}_3);\end{array}$$
 
四级四阶显式Gill公式
 
$$\{\begin{array}{l}{y}_{n+1}=y_n+\frac{h}{6}(K_1+(2-\sqrt{2})K_2+(2+\sqrt{2})K_3+K_4),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{K}_1),\\ K_3=f(x_n+\frac{1}{2}h,y_n+\frac{\sqrt{2}-1}{2}h{K}_1+(1-\frac{\sqrt{2}}{2})h{K}_2),\\ K_4=f(x_n+h,y_n-\frac{\sqrt{2}}{2}h{K}_2+(1+\frac{\sqrt{2}}{2})h{K}_3);\end{array}$$
 
示例
 
求解常微分方程:
 $$\{\begin{array}{l}\frac{dy}{dx}=x^3-\frac{y}{x},\\ y(1)=\frac{2}{5}.\end{array}$$
 要求步长为$h=0.1$，其中，精确解为
 $$y=\frac{1}{5}{x}^4+\frac{1}{5x}.$$
 
MATLAB代码
 
clear all,close all,clc
f=@(x,y)x^3-y/x;
h=0.1;
%% Euler method
x=[1:h:2];
N=size(x,2)-1;
y1=[2/5,zeros(1,N)];
for n=1:N
    y1(n+1)=y1(n)+h*f(x(n),y1(n));
end
%% Improved Euler method
y2=[2/5,zeros(1,N)];
for n=1:N
    y2(n+1)=y2(n)+h*f(x(n),y2(n));
    y2(n+1)=y2(n)+h/2*(f(x(n),y2(n))+f(x(n+1),y2(n+1)));
end
%% Standard fourth-order explicit Kutta formula
y3=[2/5,zeros(1,N)];
for n=1:N
    K1=f(x(n),y3(n));
    K2=f(x(n)+1/2*h,y3(n)+1/2*h*K1);
    K3=f(x(n)+h,y3(n)-h*K1+2*h*K2);
    y3(n+1)=y3(n)+h/6*(K1+4*K2+K3);
end
%% Three-level three-order explicit formula
y4=[2/5,zeros(1,N)];
for n=1:N
    K1=f(x(n),y4(n));
    K2=f(x(n)+1/3*h,y4(n)+1/3*h*K1);
    K3=f(x(n)+2/3*h,y4(n)+2/3*h*K2);
    y4(n+1)=y4(n)+h/4*(K1+3*K3);
end
%% Fourth-level fourth-order explicit Kutta formula
y5=[2/5,zeros(1,N)];
for n=1:N
    K1=f(x(n),y5(n));
    K2=f(x(n)+1/3*h,y5(n)+1/3*h*K1);
    K3=f(x(n)+2/3*h,y5(n)-1/3*h*K1+h*K2);
    K4=f(x(n)+h,y5(n)+h*K1-h*K2+h*K3);
    y5(n+1)=y5(n)+h/8*(K1+3*K2+3*K3+K4);
end
%% Fourth-level fourth-order explicit Gill formula
y6=[2/5,zeros(1,N)];
for n=1:N
    K1=f(x(n),y6(n));
    K2=f(x(n)+1/2*h,y6(n)+1/2*h*K1);
    K3=f(x(n)+1/2*h,y6(n)+(sqrt(2)/2-0.5)*h*K1+(1-sqrt(2)/2)*h*K2);
    K4=f(x(n)+h,y6(n)-sqrt(2)/2*h*K2+(1+sqrt(2)/2)*h*K3);
    y6(n+1)=y6(n)+h/6*(K1+(2-sqrt(2))*K2+(2+sqrt(2))*K3+K4);
end
y=1/5*x.^4+1./(5*x);  %  Exact solution

plot(x,y,'k',x,y1,'xr',x,y2,'ob','Markersize',10,'LineWidth',1.5)
axis([1 2.1 0 4.5]),xlabel x,ylabel u
legend('Exact','Euler','Improved Euler')

%plot(x,y,'k',x,y3,'*r',x,y4,'+b',x,y5,'-y',x,y6,'or','Markersize',10,'LineWidth',1.5)
%axis([1 2.1 0 6]),xlabel x,ylabel u,set(gca,'Fontsize',18)
%legend('Exact','34Kutta','33kutta','44Kutta','44Gill')

 
结果
 

科技信息SCIENCE & TECHNOLOGY INFORMATION2012 年 第 7 期 0 引言 在自然科学和经济生活的许多领域中,常常会遇到一阶常微分方程的初值问题如下：y'(t)=f(t，y)y(t0 )=y0 0 ， t0≤t (1) 这里 f 是凸区域 D 上的连续且满足 Lipschitz 条件的函数,可以证明常微分方程初值问题(1)具有唯一解. 1 常微分方程初值问题的数值解法 常微分方程初值问题的数值解法一般分为两大类，即单步法和多步法. 所谓单步法是指这类方法在计算 yn+1 时只用到前一步的 yn,这个算法的代表是 Runge-Kutta 法. Runge-Kutta 方法关于初值是稳定的, 每步的步长可以独立选取. 四阶显式 Runge-Kutta 方法是求解普通常微分方程初值问题的重要方法, 而隐式 Runge-Kutta 方法是求解刚性常微分方程初值问题的重要方法. 所谓多步法是指在计算 yn+1 时,除了用到前一步的值 yn 之外,还要用到 yn-p(p=1，2，…，k，k＞0)的值,这个算法的代表就是 Adams 方法. 1.1 龙格—库塔方法(R-K 方法) R-K 方法通过计算不同点上的函数值,并对这些函数值作线性组合,构造近似公式,利用待定系数法确定近似精度,它是非线性高阶单步法.事实上,Runge-Kutta 可以看作在(tn，tn+1)上取若干条积分曲线的若干个点的切线斜率,再进行一次(或多次)算术(或加权)平均后产生的新斜率,再按这个斜率从(tm，ym)出发,以直线带曲线向前推进一步的过 程. 由常微分方程理论,初值问题(1-1)可改写成等价的积分形式 y(t)=y0+ t t0 乙f(x，y(x))dx， 由上式及积分中值定理可得 y(tn+h)=y(tn)+hf(tn+hθ，y(tn+hθ))=yn+hk*(tn，yn，h)，0＜θ＜1. (2)但 f(tn+hθ，y(tn+hθ))未知,我们用 f 在区间［tn，tn+h］上的一些点处值的线性组合来近似它,即用下列公式代替式(2) y(tn+h)=y(tn)+h r i Σλiki, i=2，…，s, 其中 k1=fn=f(tn，yn),kj=f(tn+djh，yn+h j-1 s=1 Σβjsks)，j=2，3,…，r，dj，βjs,都是 待定的系数, 由待定系数法可以得到各阶 R-K 方法. 我们考虑 r=2 情形,应用双变量泰勒级数展式k2=f(tn+d2h，yn+β21k1h) =f(tn，yn)+h(d2 鄣 鄣t +β21k1 鄣 鄣y )f(tn，yn) +…+ h p-1 (p-1)! (d2 鄣 鄣t +β21k1 鄣 鄣y )p-1f(tn，yn)+O(hp)， 将其代入 R-K 公式,整理得 yn+1=yn+h(λ1+λ2)f(tn，yn)+h2λ2(d2ft+β21k1fy)f(tn，yn) +…+ λ2h p-1 (p-1)! (d2 鄣 鄣t +β21k1 鄣 鄣y )p-1f(tn，yn)+O(hp) 上式与 yn+1 在 yn 处的泰勒展开对比,令 h 和 h2 项的系数相同,可得 λ1+λ2=1，λ2d2= 1 2 ，λ2β21= 1 2 . d2 取不同值时可得到不同的二阶 Runge-kutta 公式. 利用上面的思想,可得经典的四阶 R-K 公式 yn+1=yn+ 1 6 h(k1+2k2+2k3+k4 ) k1=f(tn ，yn ) k2=f(tn+ 1 2 h，yn+ 1 2 hk1 ) k3=f(