28
计算方法上机报告 for mm=1:n fprintf('%9c%.f','y',mm) end for nn=1:N fprintf('\n%10.5f',x(nn)); for ll=1:n fprintf('%10.5f',y(ll,nn)); end if nn==N fprintf('\n'); end tlab 程序
附录 5 龙格—库塔法求解常微分方程初值问题的 matlab 程序
clear;clc; %======需要输入的=======% n = input('\n 请输入一阶微分方程组的个数或者高阶方程组的阶数 n：\n'); a = input('\n 请输入求解区间的下限 a:\n'); b = input('\n 请输入求解区间的上限 b:\n'); h = input('\n 请输入求解步长 h:\n'); fprintf('\n 请依次输入%.f 个方程的右端函数 fi(x,y(x)):\n',n) s = cell(n,1); for p=1:n fprintf('f%.f(x,y(x)) =',p) r = input(' ','s'); s{p} = r; end fprintf('\n 请依次输入%.f 个初值 yi(%.f):\n',n,a) y0 = zeros(n,1); for q=1:n fprintf('y%.f(%.f) =',q,a); y0(q) = input(' '); end %======需要输入的=======% N = (b-a)/h+1; x = a:h:b; K = zeros(n,4); y = zeros(n,N); y(:,1) = y0; for i=1:N-1 for a=1:n f = inline(char(s(a)),'x','y'); K(a,1) = h*f(x(i),y(:,i)); end for b=1:n f = inline(char(s(b)),'x','y'); K(b,2) = h*f(x(i)+h/2,y(:,i)+K(:,1)/2); end for c=1:n f = inline(char(s(c)),'x','y'); K(c,3) = h*f(x(i)+h/2,y(:,i)+K(:,2)/2); end for d=1:n f = inline(char(s(d)),'x','y'); K(d,4) = h*f(x(i)+h,y(:,i)+K(:,3)); end for j=1:n y(j,i+1) = y(j,i)+1/6*(K(j,1)+2*K(j,2)+2*K(j,3)+K(j,4)); end plot(x(i),y(1,i),'.k','markersize',15) hold on; end grid; fprintf('\n%10c','x')
29