【matlab】常微分方程的数值解法

实验任务

（一）常微分方程的符号计算和数值解法基本操作

1、课上例题：1、2、3、4

（二）专题实验（梯形格式）

编写梯形公式的程序，要求：

程序要有通用性，例如：

function [t,y]=tixing(f,a,b,y0,N,mytol)

% 求解形式为y’=f(t,y)的常微分方程的梯形公式

% a,b为自变量的取值范围

% y0为初始值

% N表示对区间[a,b]剖分N等份

% mytol表示允许误差, 要求默认值为1e-6，用来作为校正终止的条件

以例题4中的方程作为测试方程，画出数值解和精确解的图形，并做图例。在完成例4测试后，也可以再找其它方程作为测试方程。

（三）专题实验（改进Euler格式）

编写梯形公式的程序，要求：

（1）程序要有通用性，例如：

function [t,y]=GJ_Euler(f,a,b,y0,N)

% 求解形式为y’=f(t,y)的常微分方程的改进Euler公式

% a,b为自变量的取值范围

% y0为初始值

% N表示对区间[a,b]剖分N等份

以例4中的方程作为测试方程，画出数值解和精确解的图形，并做图例。在完成例4测试后，也可以再找其它方程作为测试方程。

（四）专题实验（变形Euler格式）

编写梯形公式的程序，要求：

（1）程序要有通用性，例如：

function [t,y]=BX_Euler(f,a,b,y0,N)

% 求解形式为y’=f(t,y)的常微分方程的变形Euler公式

% a,b为自变量的取值范围

% y0为初始值

% N表示对区间[a,b]剖分N等份

以例4中的方程作为测试方程，画出数值解和精确解的图形，并做图例。在完成例4测试后，也可以再找其它方程作为测试方程。

（五）专题实验（边值问题的有限差分方法）

编写课上的算例1、算例2，要求：

其中pfun, qfun, rfun，是方程中的系数和右端函数p(x),q(x),r(x),一般程序中不用单字母作为函数名。ceshi 是对编写的函数进行测试。myFD是有限差分方法的主程序，可以写为下列形式：

function [x,y]=myFD(alpha,beta,a,b,N)

%求常微分方程边值问题y’’+p(x)y’+q(x)y=r(x), a<x<b

y(a)=alpha, y(b)=beta

% N表示对区间[a,b]剖分N等份

%返回值x为[a,b]的剖分节点，

%返回值y为剖分节点上的数值解

（2）对课程上的算例1和算例2进行测试，每一个算例都要求画出数值解和精确解的图形，并做图例。测试的代码写在函数文件ceshi里，大家可以用其它方式编写。

（六）自选实验（可以不做）

1、编写Euler公式、向后Euler公式的程序，要求同专题实验（三）；

2、在二阶Runge-Kutta方法中选定适当的参数，编写程序，要求同专题实验（三）

实验程序和实验结果

1、实验程序：

>> y=dsolve('Dy=y^2','x')

y =

-1/(C2 + x)

2、实验程序：

>> y=dsolve('x*D2y-3*Dy=x^2','y(1)=0','y(5)=0','x')

y =

(31*x^4)/468 - x^3/3 + 125/468

3、实验程序：

>> y=dsolve('Df=f+g','Dg=-f+g','f(0)=1','g(0)=2','t')

y =

包含以下字段的 struct:

g: [1×1 sym]

f: [1×1 sym]

>> y.f

ans =

exp(t)*cos(t) + 2*exp(t)*sin(t)

>> y.g

ans =

2*exp(t)*cos(t) - exp(t)*sin(t)

4、实验程序：

[x,y]=ode23(@funst,[0,1],pi/2);

yy=(x+pi/2)./cos(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','解析解')

function yp=funst(x,y)

yp=y*tan(x)+sec(x);

end

（二）专题实验（梯形格式）

1、实验程序：

function [x,y]=tixing(f,a,b,y0,N,mytol)

% 求解形式为y’=f(t,y)的常微分方程的梯形公式

% a,b为自变量的取值范围

% y0为初始值

% N表示对区间[a,b]剖分N等份

% mytol表示允许误差, 要求默认值为1e-6，用来作为校正终止的条件

%nargin表示参数个数

if nargin<=5

mytol=1e-6;

end

h=(b-a)/N;

x=a:h:b;

y(1)=y0;

for i=1:length(x)-1

y(i+1)=y(i)+h*f(x(i),y(i));

y1=y(i+1)+1;

while abs(y(i+1)-y1)>=mytol

y1=y(i+1);

y(i+1)=y(i)+h/2*(f(x(i),y(i))+f(x(i+1),y1));

end

实验程序：

先建立函数文件funst.m

function yp=funst(x,y)

yp=y*tan(x)+sec(x);

end

求解程序：

[x,y]=tixing(@funst,0,1,pi/2,10);

yy=(x+pi/2)./cos(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','解析解')

（三）专题实验（改进Euler格式）

1、实验程序：

function [x,y]=GJ_Euler(f,a,b,y0,N)

% 求解形式为y’=f(t,y)的常微分方程的梯形公式

% a,b为自变量的取值范围

% y0为初始值

% N表示对区间[a,b]剖分N等份

% mytol表示允许误差, 要求默认值为1e-6，用来作为校正终止的条件

%nargin表示参数个数

if nargin<=5

mytol=1e-6;

else

error

end

h=(b-a)/N;

x=a:h:b;

y(1)=y0;

for i=1:length(x)-1

y(i+1)=y(i)+h*f(x(i),y(i));

y(i+1)=y(i)+h/2*(f(x(i),y(i))+f(x(i+1),y(i+1)));

end

2、实验程序：

先建立函数文件funst.m

function yp=funst(x,y)

yp=y*tan(x)+sec(x);

end

求解程序:

[x,y]=GJ_Euler(@funst,0,1,pi/2,10);

yy=(x+pi/2)./cos(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','解析解')

专题实验（变形Euler格式）

实验程序：

function [x,y]=BX_Euler(f,a,b,y0,N)

% 求解形式为y’=f(t,y)的常微分方程的梯形公式

% a,b为自变量的取值范围

% y0为初始值

% N表示对区间[a,b]剖分N等份

% mytol表示允许误差, 要求默认值为1e-6，用来作为校正终止的条件

%nargin表示参数个数

if nargin<=5

mytol=1e-6;

else

error

end

h=(b-a)/N;

x=a:h:b;

y(1)=y0;

for i=1:length(x)-1

y(i+1)=y(i)+h*f(x(i)+h/2,y(i)+h/2*f(x(i),y(i)));

end

实验程序：

先建立函数文件funst.m

function yp=funst(x,y)

yp=y*tan(x)+sec(x);

end

求解程序：

[x,y]=BX_Euler(@funst,0,1,pi/2,10);

yy=(x+pi/2)./cos(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','解析解')

（五）专题实验（边值问题的有限差分方法）

1、实验程序：

function [x,y]=myFD(alpha,beta,a,b,N)

%求常微分方程边值问题y''+p(x)y’+q(x)y=r(x),a<x<b,y(a)=alpha,y(b)=beta

% N表示对区间[a,b]剖分N等份

%返回值x为[a,b]的剖分节点

%返回值y为剖分节点上的数值解

h=(b-a)/N;

x=a:h:b;

y=zeros(length(x),1);

y(1)=alpha; y(length(x))=beta;

f=zeros(N-1,1);

A=zeros(N-1,N-1);

for i=1:N-1

if i==1

f(i,1)=h^2*rfun(x(2))-(1-h/2*pfun(x(2)))*alpha;

A(1,1)=-2+h^2*qfun(x(2));

A(1,2)=1+h/2*pfun(x(2));

elseif i==N-1

f(i,1)=h^2*rfun(x(N))-(1+h/2*pfun(x(N)))*beta;

A(N-1,N-2)=1-h/2*pfun(x(N));

A(N-1,N-1)=-2+h^2*qfun(x(N));

else

f(i,1)=h^2*rfun(x(i+1));

A(i,i-1)=1-h/2*pfun(x(i+1));

A(i,i)=-2+h^2*qfun(x(i+1));

A(i,i+1)=1+h/2*pfun(x(i+1));

end

y(2:N,1)=A\f;

end

2、实验程序：

function px=pfun(x)

px=3;

end

function qx=qfun(x)

qx=2;

end

function rx=rfun(x)

rx=sin(x)+3*cos(x);

end

[x,y]=myFD(0,1,0,pi/2,10);

yy=sin(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','精确解')

实验程序：

function px=pfun(x)

px=2+x^2;

end

function qx=qfun(x)

qx=1+2*x^2;

end

function rx=rfun(x)

rx=-sin(x)+(2+x^2)*cos(x)+(1+2*x^2)*sin(x);

end

[x,y]=myFD(0,1,0,pi/2,10);

yy=sin(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','精确解')

【matlab】常微分方程的数值解法相关推荐

二阶常微分方程的数值解法（中心差分法和有限体积法）
二阶常微分方程的数值解法(中心差分法和有限体积法) 这里我们介绍中心差分法和有限体积法求解方程. 题目: 用差分法的中心差分格式和有限体积法求解两点边值问题 u′′−α(2x−1)u′−2αu=0,0 ...
常微分方程数值解的c语言程序,常微分方程的数值解法一阶常微分方程数值解的C语言编程实现.doc...
常微分方程的数值解法一阶常微分方程数值解的C语言编程实现导读:就爱阅读网友为您分享以下"一阶常微分方程数值解的C语言编程实现"资讯,希望对您有所帮助,感谢您对92的支持! 一阶 ...
【数理知识】《数值分析》李庆扬老师-第9章-常微分方程初值问题数值解法
第8章回到目录无第9章-常微分方程初值问题数值解法 9.1 引言利普希茨 (Lipschitz) 条件 / 利普希茨常数定理1 解的存在唯一性定理定理2 解对初值依赖的敏感性 9.2 简单 ...
[常微分方程的数值解法系列五] 龙格-库塔(RK4)法
龙格-库塔法简介基本思想具体方法一阶二阶求解参数特殊二阶三阶高阶步长选择例子在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其 ...
一阶常微分方程的数值解法（二阶显式、隐式 Adams 公式及 Milne 方法）
一阶常微分方程的数值解法这里我们介绍二阶显式.隐式 Adams 公式及 Milne 方法求解方程. 题目: 对初值问题 u ′ = u − t 2 , 0 ≤ t ≤ 1 , u ( 0 ) = 0 ...
常微分方程初值问题数值解法[完整公式]（Python）
目录 1.概述 (1)常微分方程初值问题数值解法 (2)解题步骤 (3)数值微分解法 (4)数值积分解法 2.所有知识点代码 3.结果---以三阶Runge-Kutta公式为例(其他的类似) 1.概述 ...
[常微分方程的数值解法系列三] 改进欧拉法(预估校正法)
改进欧拉法简介预估-校正截断误差例子在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其中主要就是求解微分方程.但之前求解微分方程的解析方法主 ...
[常微分方程的数值解法系列二] 欧拉法
欧拉法简介几何意义证明泰勒展开近似求导近似积分近似几种欧拉方式向前欧拉公式向后欧拉公式梯形公式中点公式截断误差求解过程向前欧拉公式例子向前欧拉公式在惯性导航以及VIO ...
[常微分方程的数值解法系列四] 中值法
中值法简介具体步骤截断误差例子在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其中主要就是求解微分方程.但之前求解微分方程的解析方法主要是应 ...
二阶龙格库塔公式推导_[常微分方程的数值解法系列五] 龙格-库塔(RK4)法
在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其中主要就是求解微分方程.但之前求解微分方程的解析方法主要是应用于一些简单和特殊的微分方程求解中,对于 ...

【matlab】常微分方程的数值解法

【matlab】常微分方程的数值解法相关推荐

最新文章

热门文章