matlab数值分析程序--高等数学,数值代数的matlab实现-文字版, matlab电子书, 和matlab 有关的电子书:

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

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

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

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]上的数值解','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,:);fxy=fxy(1:k,:);

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

wch=wch(1:k,:);

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

用常用四阶龙格-库塔方法求解常微分方程初值问题(10.5)的数值解的MATLAB主程序2

function [x,y]=RKc4(x0,b,y0,h)

for x=x0:h:b

    k1=funfcn(x0,y0);

k2=funfcn(x0+h/2,y0+h*k1/2);

    k3=funfcn(x0+h/2,y0+h*k2/2);

k4=funfcn(x0+h,y0+h*k3);

    y=y0+h*(k1+2*k2+2*k3+k4)/6;

x=x0+h;[x,y], x0=x;y0=y;

end

[x,y];

 

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

的数值解,h=1/4,并计算与精确解的误差,在同一图形窗口画出例10.5.1、例10.5.2和例10.5.3精确解和数值解的图形.

解   方法1在MATLAB工作窗口输入下面的程序

>> x0=0;b=2; y0=0;h=1/4;

subplot(3,1,1)

C=[1/4,3/4,2/3,2/3];

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

subplot(3,1,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];

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

subplot(3,1,3)

c1=1/6;c2=2/6;c3=2/6;c4=1/6;a2=1/2; a3=1/2; 

a4=1;b21=1/2;b31=0;b32=1/2;

b41=0;b42=0;b43=1;

C=[c1,c2,c3, c4,a2, a3, a4,b21,b31,b32,b41,b42,b43];

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

将运行后计算的结果(略),画出精确解和数值解的图形.

 

方法2  在MATLAB工作窗口输入下面的程序

>> x0=0;b=2; y0=0;h=1/4;

[x,y]=RKc4(x0,b,y0,h)

运行后输出结果(略).

欢迎转载,转载请注明来自一手册:http://yishouce.com/book/3/3100504.html
友情链接It题库(ittiku.com)| 版权归yishouce.com所有| 友链等可联系 admin#yishouce.com|粤ICP备16001685号-1