[原创]龙格-库塔法求解微分方程

wangzi629

function yy=weifen_lgkt(A,B,x0,y0,fx)
% yy=weifen_lgkt(A,B,x0,y0,fx)
% 系统状态方程
% dy/dx-Ay=B*fx 微分方程,A,B为参数
% fx x的函数,可选参数,如不输入该参数, 则要在该函数体内编辑该函数
% x0,y0 状态变量初始值%注：x0,y0可以为列向量,求解微分方程组
N=1000; %计算步数
h=0.002; %计算步长
x=x0; y=y0;
%---------------------
for k=1:N
%求解微分方程
%fx(:,k)=x0+(k-1)*h; %例子1，在此选f(x)=x
k0=A*y+B*fx(:,k); %四阶龙格-库塔法则求解微分方程
k1=A*(y+h*k0/2)+B*fx(:,k); k2=A*(y+h*k1/2)+B*fx(:,k);
k3=A*(y+h*k2)+B*fx(:,k); y=y+(k0+2*k1+2*k2+k3)*h/6;
%--------------------------------------------------------------------------------------%
yy(k)=y; %输出
xx(k)=x0+k*h;

%y_jiexi(k)=2*fx(:,k)-2+4*exp(-fx(:,k)); %例子1的解析解
%y_e(k)=y_jiexi(k)-yy(k); %例子1解析解与数值解的误差
end
plot([x0 xx],[y0 yy])

*****************************************************************************************
程序完毕。
说明：如求解方程dy/dx+y=2x 初值x0=0 y0=2 则选A=-1 B=2 x0=0 y0=2
fx=x 见程序内例1实现方法 ，当然也可把fx当作参数输进去。
调用函数即可（例1需把%fx(:,k)=x0+(k-1)*h;前的%去掉）
A=-1; B=2; x0=0; y0=2;
yy=weifen_lgkt(A,B,x0,y0)

就可以求得该微分方程

[ 本帖最后由 ChaChing 于 2010-6-29 08:26 编辑 ]

hunter_009

请问，楼主，你的方法与ode23或者ode45有什么不同呢？又各有什么优缺点呢？

fandalei

请问楼主积分步长如何选择
楼主好象不怎么来哦,那请happy 老师来给我们讲讲了
楼主,A,B, C是矩阵,怎么办?

[ 本帖最后由 ChaChing 于 2010-6-29 08:41 编辑 ]

wangzi629

步长h对数值方法的稳定性影响很大。由于内容涉及太多，可以找本《数值计算方法》或《数值分析》方面的书看看。上面有关于微分方程求解方面的知识。比如dy/dx=ry.当r<0时，四阶龙格-库塔的绝对稳定域为-2.78<rh<0.即h<2.78/(-r)时此方法才稳定
A，B，C是矩阵时可以照常使用，但要注意各矩阵维数应满足矩阵乘法。

[ 本帖最后由 ChaChing 于 2010-6-29 08:29 编辑 ]

fandalei

A,B, C是矩阵,fx(:,k)怎么办,它是随时间变化的函数或者是不变的;
楼主说的步长是1阶微分方程,如果是2阶微分方程怎么办?

fx(:,k)怎么办,它是随时间变化的函数

斑竹，我的帖子可以设为精华吗？

[ 本帖最后由 ChaChing 于 2010-6-29 08:35 编辑 ]

grta

有没有变步长的程序啊

odovo

ode45就是变步长的

wenyue

A,B,C是矩阵时请高手赐教

AaronSpark

原帖由 wenyue 于 2006-9-21 23:33 发表
A,B,C是矩阵时请高手赐教

换成矩阵运算就可以了
建议到matlab赏析版找找happy发的程序
和ode45是比较完善的，可以算你的这种情况

原帖由 fandalei 于 2006-6-27 08:39 发表
<P>请问楼主积分步长如何选择</P>

找本数值分析的书看看吧

[ 本帖最后由 ChaChing 于 2010-6-29 08:38 编辑 ]

AaronSpark

原帖由 hunter_009 于 2006-6-28 22:02 发表
请问，楼主，你的方法与ode23或者ode45有什么不同呢？又各有什么优缺点呢？

阶数不一样

ode45精度高，计算慢一些
ode23精度低，计算速度快
一般采用ode45比较多

原帖由 fandalei 于 2006-7-1 15:19 发表
fx(:,k)怎么办,它是随时间变化的函数

那就加上时间表达式

[ 本帖最后由 ChaChing 于 2010-6-29 08:38 编辑 ]

ak47kill

其中的fx是怎么编辑的啊
如果fx=exp(x)，应该怎么办啊

FtpAdmin

原帖由 hunter_009 于 2006-6-28 22:02 发表
请问，楼主，你的方法与ode23或者ode45有什么不同呢？又各有什么优缺点呢？

ode23和ode45都是变步长的
楼主的程序是定步长的

FtpAdmin

原帖由 ak47kill 于 2006-10-22 16:11 发表
其中的fx是怎么编辑的啊
如果fx=exp(x)，应该怎么办啊

写一个这个函数的子程序就行了

FtpAdmin

原帖由 fandalei 于 2006-6-29 19:07 发表
<P>A,B, C是矩阵,fx(:,k)怎么办,它是随时间变化的函数或者是不变的;<BR>楼主说的步长是1阶微分方程,如果是2阶微分方程怎么办?</P>

楼主的程序缺乏通用性，建议采用本版happy提供的几个定步长RK法程序

yanyongju

ode23是解决刚性问题的，而ode45是主要解决非刚性问题的，在运算时间是ode45要常于ode23，ode45是四五阶的，精度要高于ode23，如果系数是矩阵，则涉及到了解藕，需要将矩阵化为对角阵。