第十章 常微分方程(组)求解(6)

第三篇 第十章 常微分方程（组）求解

h?y?y?(k1?3k2),n?n?14? (10.39) ?k1?f(xn,yn),?22?k2?f(xn?h,yn?hk1).33?可以证明，在[xn,xn?1]内只取2点的龙格-库塔公式精度最高为2阶。

用二阶龙格—库塔方法求解常微分方程初值问题（10.5）的数值解的MATLAB主程序 输入量：funfcn是函数f(x,y)，fun是初值问题的精确解函数，x0和y0是初值问题y(x0)?y0,b是自变量x的最大值，C?[c1,c2,a2,b21]是(10.36)式的参数组成的数组，h是步长。

输出量：向量X的元素是由自变量x的取值组成，数组Y的元素是利用(10.33)式、(10.36)式求出常微分方程初值问题(10.5)在X的元素处的数值解，fxy是在向量X的元素处的精确解，n是自变量x取值的序号，wucha(k+1)=norm(Y(k+1)-Y(k)),wch(k)=norm(fxy(k)-Y(k)),并画出数值解向量Y的图形。

function

[k,X,Y,fxy,wch,wucha,P]=RK2(funfcn,fun,x0,b,C,y0,h)

x=x0; y=y0;p=128; n=fix((b-x0)/h); fxy=zeros(p,1); wucha=zeros(p,1); wch=zeros(p,1); X=zeros(p,1);

Y=zeros(p,length(y)); k=1; X(k)=x; Y(k,:)=y'; % 绘图.

clc,x,h,y %计算

%fxy=fxy(:); for k=2:n+1

x=x+h;a2=C(3);b21=C(4); c1=C(1); c2=C(2); x1=x+a2*h;

k1=feval(funfcn,x,y); y1=y+b21*h*k1;

k2=feval(funfcn,x1,y1); fxy(k)=feval(fun,x); y=y+h*(c1*k1+c2*k2); X(k)=x;

Y(k,:)=y;k=k+1;

plot(X,Y,'mh',X,fxy,'bo')

grid,xlabel('自变量 X'),ylabel('因变量 Y')

legend('用二阶龙格-库塔方法计算dy/dx=f(x,y)，y(x0)=y0在

[x0,b]上的数值解','y/dx=f(x,y)，y(x0)=y0的精确解y=f(x)')

end

%计算误差. for k=2:n+1

wucha(k)=norm(Y(k-1)-Y(k)); wch(k)=norm(fxy(k)-Y(k)); end

195．

第三篇 第十章 常微分方程（组）求解

X=X(1:k); Y=Y(1:k,:);

fxy=fxy(1:k,:); n=1:k;

wucha=wucha(1:k,:); wch=wch(1:k,:);

P=[n',X,Y,fxy,wch,wucha];

例10.5.1 用二阶龙格—库塔方法求初值问题

2xy?dy?1?,0?x?2, ?2dx1?x???yx?0?0的数值解，取c1?1/4,c2?3/4,a2?b21?2/3，h?1/4，并计算与精确解的误差，画出精确解和数值解的图形.

解 （1）先求精确解。输入程序：

>>y=dsolve('(Dy)+(2*x*y)/(1+x^2)-1=0','y(0)=0','x') (2) 编写并保存名为funfcn.m和fun.m的文件如下 function f=funfcn(x,y) f=1-(2*x*y)/(1+x^2); function y=fun(x)

y=(x+1/3*x^3)/(1+x^2);

(3)在MATLAB工作窗口输入下面的程序

>> x0=0;b=2;C=[1/4,3/4,2/3,2/3];y0=0;h=1/4;

[k,X,Y,fxy,wch,wucha,P]=RK2(@funfcn,@fun,x0,b,C,y0,h) ④将运行后计算的结果列入（略），画出精确解和数值解的图形.

10.5.3 三阶龙格-库塔方法及其MATLAB程序 （一） 三阶龙格-库塔方法 要进一步提高精度，必须取更多的点，如果在区间[xn,xn?1]内取r?3个点，

由(10.32)式和(10.31)式，得三阶龙格-库塔方法的形式为

?yn?1?yn?h(c1k1?c2k2?c3k3),?k?f(x,y),nn?1 (10.40) ?k?f(x?ah,y?bhk),n2n211?2?k3?f?xn?a3h,yn?h(b31k1?b32k2)?,??c1?c2?c3?1,??a2c2?a3c3?1,2??212c3?,?a2c2?a3其中?3 (10.41)

?1?a2b32c3?,6??a2?b21,??a3?b1?b32.因为方程组(10.41)含八个未知数和六个方程，所以有多个不同的解，因此有多种不同的三阶龙格-库塔公式。取

1211c1?,c2?,c3?,a2?b21?,a3?1,b31??1,b32?2,，则(10.40)式为常用

6362 196．

第三篇 第十章 常微分方程（组）求解

的三阶龙格-库塔公式。

用三阶龙格-库塔方法求解常微分方程初值问题（10.5）的数值解的MATLAB主程序 输入量：funfcn是函数f(x,y)，fun是初值问题的精确解函数，x0和y0是初值问题y(x0)?y0,b是自变量x的最大值，C?[c1,c2,c3,a2,a3,b21,b31,b32]是(10.40)式的参数组成的数组，h是步长。

输出量：向量X的元素是由自变量x的取值组成，数组Y的元素是利用(10.40)式、(10.41)式求出常微分方程初值问题(10.5)在X的元素处的数值解，fxy是在向量X的元素处的精确解，n是自变量x取值的序号，wucha(k+1)=norm(Y(k+1)-Y(k)),wch(k)=norm(fxy(k)-Y(k)),并画出数值解向量Y的图形。

function

[k,X,Y,fxy,wch,wucha,P]=RK3(funfcn,fun,x0,b,C,y0,h)

x=x0; y=y0;p=128; n=fix((b-x0)/h); fxy=zeros(p,1); wucha=zeros(p,1); wch=zeros(p,1); X=zeros(p,1);

Y=zeros(p,length(y)); k=1; X(k)=x; Y(k,:)=y'; % 绘图.

clc,x,h,y %计算

%fxy=fxy(:); for k=2:n+1

x=x+h;a2=C(4); a3=C(5);b21=C(6); b31=C(7);

b32=C(8);c1=C(1); c2=C(2); c3=C(3); x1=x+a2*h; x2=x+a3*h;

k1=feval(funfcn,x,y); y1=y+b21*h*k1;

k2=feval(funfcn,x1,y1); y2=y+b31*h*k1+b32*h*k2; k3=feval(funfcn,x2,y2); fxy(k)=feval(fun,x);

y=y+h*(c1*k1+c2*k2+c3*k3); X(k)=x; Y(k,:)=y; k=k+1;

plot(X,Y,'rp',X,fxy,'bo'),

grid,xlabel('自变量 X'), ylabel('因变量 Y')

legend('用三阶龙格-库塔方法计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]

上的数值解','y/dx=f(x,y)，y(x0)=y0的精确解y=f(x)')

end

%计算误差. for k=2:n+1

wucha(k)=norm(Y(k-1)-Y(k)); wch(k)=norm(fxy(k)-Y(k)); end

X=X(1:k); Y=Y(1:k,:);

197．

第三篇 第十章 常微分方程（组）求解

fxy=fxy(1:k,:); n=1:k;

wucha=wucha(1:k,:); wch=wch(1:k,:);

P=[n',X,Y,fxy,wch,wucha];

例10.5.2 用常用的三阶龙格-库塔公式求初值问题

2xy?dy,0?x?2,??1? 1?x2?dx??yx?0?0的数值解，取c1?1/6,c2?4/6,c3?1/6, a2?b21?1/2，h=1/4,a3?1,b31??1,b32?2，

并计算与精确解的误差，画出精确解和数值解的图形.

解 (1) 编写并保存名为funfcn.m和fun.m的文件如下 function f=funfcn(x,y) f=1-(2*x*y)/(1+x^2); function y=fun(x)

y=(x+1/3*x^3)/(1+x^2);

（2）在MATLAB工作窗口输入下面的程序

>> x0=0;b=2; c1=1/6; c2=4/6;c3=1/6;a2=1/2;

a3=1;b21=1/2;b31=-1;b32=2;

C=[c1,c2,c3,a2, a3,b21,b31,b32]; y0=0;h=1/4;

[k,X,Y,fxy,wch,wucha,P]=RK3(@funfcn,@fun,x0,b,C,y0,h)

将运行后计算的结果（略），画出精确解和数值解的图形.

10.5.4 四阶龙格-库塔方法及其MATLAB程序 （一） 四阶龙格-库塔方法

如果在区间[xn,xn?1]内取r?4个点，由(10.32)式和(10.31)式，得四阶龙格-库塔方法的形式为

?y?y?h(ck?ck?ck?ck),n11223344?n?1?k1?f(xn,yn),? (10.42) ?k2?f(xn?a2h,yn?b21hk1),?k?fx?ah,y?h(bk?bk),?n3n311322??3?k?f?x?ah,y?h(bk?bk?bk)?,n4n411422433?4其中待定系数ci,ai,bij公13个，经过与推导二阶龙格-库塔公式类似，但更复

杂的计算，得到使局部误差y(xn?1)?yn?1?O(h5)的11个方程。取即满足这些

方程、又较简单的一组ci,ai,bij,可得

198．

第三篇 第十章 常微分方程（组）求解

?h?yn?1?yn?(k1?2k2?2k3?k4),6??k1?f(xn,yn),?h1? (10.43) ?k2?f(xn?,yn?hk1),22??k2?h?k?fx?,y?h,nn?3??22?????k4?f?xn?h,yn?hk3?.这就是常用的四阶龙格-库塔公式。

用一般四阶龙格-库塔方法（10.42）式求解常微分方程初值问题（10.5）的数值解的MATLAB主程序1

输入量：funfcn是函数f(x,y)，fun是初值问题的精确解函数，x0和y0是初值问

是自变量的最大值，y(x0)?y0,bxC?[c1,c2,c3,c4,a2,a3,a4,b21,b31,b32,b41,b42,b43]是(10.42)式的参数组成的数组，h是步长。 输出量：向量X的元素是由自变量x的取值组成，数组Y的元素是利用(10.42)式求出常微分方程初值问题(10.5)在X的元素处的数值解，fxy是在向量X的元素处的精确解，n是自变量x取值的序号，

wucha(k+1)=norm(Y(k+1)-Y(k)),wch(k)=norm(fxy(k)-Y(k)),并画出数值解向量Y的图题

形。

function [k,X,Y,fxy,wch,wucha,P]=RK4

(funfcn,fun,x0,b,C,y0,h)

x=x0; y=y0;p=128; n=fix((b-x0)/h); fxy=zeros(p,1); wucha=zeros(p,1); wch=zeros(p,1); X=zeros(p,1); Y=zeros(p,length(y)); k=1; X(k)=x; Y(k,:)=y'; % 绘图.

clc,x,h,y %计算

%fxy=fxy(:); for k=2:n+1

x=x+h;a2=C(5); a3=C(6);

a4=C(7);b21=C(8); b31=C(9); b32=C(10); b41=C(11);

b42=C(12); b43=C(13);c1=C(1); c2=C(2); c3=C(3); c4=C(4); x1=x+a2*h; x2=x+a3*h; x3=x+a4*h;

k1=feval(funfcn,x,y); y1=y+b21*h*k1;

k2=feval(funfcn,x1,y1); y2=y+b31*h*k1+b32*h*k2; k3=feval(funfcn,x2,y2);

y3=y+b41*h*k1+b42*h*k2+b43*h*k3; k4=feval(funfcn,x3,y3); fxy(k)=feval(fun,x);

y=y+h*(c1*k1+c2*k2+c3*k3+c4*k4); X(k)=x; Y(k,:)=y;

k=k+1; plot(X,Y,'rp',X,fxy,'bo'), grid xlabel('自变量 X'), ylabel('因变量 Y')

legend('用四阶龙格-库塔方法计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]

199．

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库第十章 常微分方程(组)求解(6)在线全文阅读。