matlab编的4阶龙格库塔法解微分方程的程序.docVIP

matlab编的4阶龙格库塔法解微分方程的程序 2010-03-10 20:16 function varargout=saxplaxliu(varargin)clc,clearx0=0;xn=1.2;y0=1;h=0.1;[y,x]=lgkt4j(x0,xn,y0,h); n=length(x);fprintf( i?? x(i)?? y(i)\n);for i=1:n??? fprintf(%2d %4.4f %4.4f\n,i,x(i),y(i));endfunction z=f(x,y)?? z=-2*x*y^2;function [y,x]=lgkt4j(x0,xn,y0,h) x=x0:h:xn;n=length(x);y1=x;y1(1)=y0;for i=1:n-1??? K1=f(x(i),y1(i));??? K2=f(x(i)+h/2,y1(i)+h/2*K1);??? K3=f(x(i)+h/2,y1(i)+h/2*K2);??? K4=f(x(i)+h,y1(i)+h*K3);??? y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);end?? y=y1;?? 结果： i?? x(i)?? y(i)1 0.0000 1.00002 0.1000 0.99013 0.2000 0.96154 0.3000 0.91745 0.4000 0.86216 0.5000 0.80007 0.6000 0.73538 0.7000 0.67119 0.8000 0.609810 0.9000 0.552511 1.0000 0.500012 1.1000 0.452513 1.2000 0.4098 龙格库塔法 一、基本原理：龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。对于一阶精度的欧拉公式有： yi+1=yi+h*K1 K1=f(xi,yi) 当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进欧拉公式： yi+1=yi+h*( K1+ K2)/2 K1=f(xi,yi) K2=f(xi+h,yi+h*K1) 依次类推，如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法： yi+1=yi+h*( K1+ 2*K2 +2*K3+ K4)/6 K1=f(xi,yi) K2=f(xi+h/2,yi+h*K1/2) K3=f(xi+h/2,yi+h*K2/2) K4=f(xi+h,yi+h*K3) 通常所说的龙格-库塔法是指四阶而言的，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式 （1） 计算公式（1）的局部截断误差是 。 龙格-库塔法具有精度高，收敛，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算 在一些点上的值，如四阶龙格-库塔法每计算一步需要计算四次 的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”。 二、小程序 #includestdio.h #includemath.h #define f(x,y) (-1*(x)*(y)*(y)) void main(void) { double a,b,x0,y0,k1,k2,k3,k4,h; int n,i; printf(input a,b,x0,y0,n:); scanf(%lf%lf%lf%lf%d,a,b,x0,y0,n); printf(x0\ty0\tk1\tk2\tk3\tk4\n); for(h=(b-a)/n,i=0;i!=n;i++) { k1=f(x0,y0); k2=f(x0+h/2,y0+k1*h/2); k3=f(x0+h/2,y0+k2*h/2); k4=f(x0+h,y0+h*k3); printf(%lf\t%lf\t,x0,y0); printf(%lf\t%lf\t,k1,k2); printf(%lf\t%lf\n,k3,k4); y0+=h*(k1+2*k2+2*k3+k4)/6; x0+=h; } printf(xn=%lf\tyn=%lf\n,x0,y0); } 运行结果： input a,b,x0,y0,n:0 5 0 2 20 x0? ?? ?y0? ?? ?k