微分方程組的龍格庫塔公式求解matlab版

微分方程組的龍格-庫塔公式求解matlab版

南京大學(xué) 王尋

1. 一階常微分方程組

考慮方程組

? ? ? ?

?y'?f x,y,z , y x ?y

? 0 0

? ? ? ?

z'?g x,y,z , z x ?z

? 0 0

其經(jīng)典四階龍格-庫塔格式如下：

對于n=0,1,2,...,計(jì)算

? h

y ?y ? K ?2K ?2K ?K? ?

? n?1 n 6 1 2 3 4

? h

?z ?z ? L ?2L ?2L ?L? ?

n?1 n 1 2 3 4

? 6

其中

?K ?f x ,y ,z , L ?g x ,y ,z? ? ? ?

1 n n n 1 n n n

? ? h hK1 hL1? ? h hK1 hL1?

?K ?f x ? ,y ? ,z ?2 ? n n n ?,L ?g x ? ,y ? ,z ?2 ? n n n ?

? ? 2 2 2 ? ? 2 2 2 ?

?

?K ?f x ? ,y ?3 ?? n h n hK2 ,z ?n hL2 ??,L ?g x ? ,y ?3 ?? n h n hK2 ,z ?n hL2 ??

? ? 2 2 2 ? ? 2 2 2 ?

?K ?f x ?h,y ?hK ,z ?hL ,L ?g x ?h,y ?hK ,z ?hL? ? ? ?

? 4 n n 3 n 3 4 n n 3 n 3

下面給出經(jīng)典四階龍格-庫塔格式解常微分方程組 matlab通用程序：

%marunge4s.m

%用途：4階經(jīng)典龍格庫塔格式解常微分方程組y'=f(x,y),y(x0)=y0

%格式：[x,y]=marunge4s(dyfun,xspan,y0,h)

%dyfun為向量函數(shù)f(x,y),xspan為求解區(qū)間[x0,xn],

%y0為初值向量，h為步長，x返回節(jié)點(diǎn)，y返回?cái)?shù)值解向量

function [x,y]=marunge4s(dyfun,xspan,y0,h)

x=xspan(1):h:xspan(2);

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

y(:,1)=y0(:);

for n=1:(length(x)-1)

k1=feval(dyfun,x(n),y(:,n));

k2=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k1);

k3=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k2);

k4=feval(dyfun,x(n+1),y(:,n)+h*k3);

y(:,n+1)=y(:,n)+(h/6).*(k1+2*k2+3*k3+k4);

end

如下為例題：

例1：取h=0.02，利用程序marunge4s.m求剛性微分方程組

? ?

?y'??0.01y?99.99z, y 0 ?2,

?

? ?

?z'??100z, z 0 ?1

的數(shù)值解，其解析解為： ?0.01x ?100x ?100x

y ?e ?e ,z?e

M dyfun.m

解：首先編寫 函數(shù)

%dyfun.m

function f=dyfun(t,y)

f(1)=-0.01*y(1)-99.99*y(2);

f(2)=-100*y(2);

f=f(:);

與50位技術(shù)專家面對面20年技術(shù)見證，附贈技術(shù)全景圖

總結(jié)

以上是生活随笔為你收集整理的龙格库塔法matlab求解微分方程组,微分方程组的龙格库塔公式求解matlab版.pdf的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。