龙格库塔法的编程

龙格库塔法的编程龙格库塔法的编程 #include #include /*n表示几等分，n+1表示他输出的个数*/ int RungeKutta(double y0,double a,double b,int n,double *x,double *y,int style,double (*function)(double，double)) { double h=(b-a)/n,k1,k2,k3,k4; int i; // x=(double*)malloc((n+1)*sizeof(double)); // y...

龙格库塔法的编程 #include #include /*n 表示几等分，n+1表示他输出的个数*/ int RungeKutta(double y0,double a,double b,int n,double *x,double *y,int style,double (*function)(double，double)) { double h=(b-a)/n,k1,k2,k3,k4; int i; // x=(double*)malloc((n+1)*sizeof(double)); // y=(double*)malloc((n+1)*sizeof(double)); x[0]=a; y[0]=y0; switch(style) { case 2: for(i=0;i #include /*n表示几等分，n+1表示他输出的个数*/ int RungeKutta(double y0,double a,double b,int n,double *x,double *y,int style,double (*function)(double,double)) { double h=(b-a)/n,k1,k2,k3,k4; int i; // x=(double*)malloc((n+1)*sizeof(double)); // y=(double*)malloc((n+1)*sizeof(double)); x[0]=a; y[0]=y0; switch(style) { case 2: for(i=0;i 工程中很多的地方用到龙格库塔求解微分方程的数值解，龙格库塔是很重要的一种方法，尤其是四阶的，精确度相当的高。此代码只是演示求一个微分方程：的解，要求解其它的微分方程，可以自己定义借口函数，退换程序里面的函数： float f(float , float)即可； CODE: /* /*编程者：刘艮平 /*完成日期：2004.5.29 /*e.g:y'=y-2x/y,x∈[0,0.6] /*　　y(0)=1 /*使用经典四阶龙格-库塔算法进行高精度求解　 /* y(i＋1)＝yi＋( K1+ 2*K2 +2*K3+ K4)/6 /* K1=h*f(xi,yi) /* K2=h*f(xi+h/2,yi+K1/2) /* K3=h*f(xi+h/2,yi+K2/2) /* K4=h*f(xi+h,yi+K3) */ #include #include float f(float den,float p0) //要求解的微分方程的右部的函数 e.g: y-2x/y { float rus; // den=w/(W0+sl); // rus=k*A*f/Wp*sqrt(RTp)-(k-1)*..... rus=p0-2*den/p0; return(rus); } //float fden() //{ //} void main() { float x0; //范围上限 int x1; //范围下限 float h; //步长 int n; //计算出的点的个数 float k1,k2,k3,k4; float y[3]; //用于存放计算出的常微分方程数值解 int i=0; int j; //以下为函数的接口 printf("intput x0:"); scanf("%f",&x0); printf("input x1:"); scanf("%f",&x1); printf("input y[0]:"); scanf("%f",&y[0]); printf("input h:"); scanf("%f",&h); // 以下为核心程序 n=((x1-x0)/h); n=3; for(j=0;j 　　#include 　　#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 k1 k2 k3 k4 　　0.000000 2.000000 -0.000000 -0.500000 -0.469238 　　-0.886131 　　0.250000 1.882308 -0.885771 -1.176945 -1.129082 　　-1.280060 　　0.500000 1.599896 -1.279834 -1.295851 -1.292250 　　-1.222728 　　0.750000 1.279948 -1.228700 -1.110102 -1.139515 　　-0.990162 　　1.000000 1.000027 -1.000054 -0.861368 -0.895837 　　-0.752852 　　1.250000 0.780556 -0.761584 -0.645858 -0.673410 　　-0.562189 　　1.500000 0.615459 -0.568185 -0.481668 -0.500993 　　-0.420537 　　1.750000 0.492374 -0.424257 -0.361915 -0.374868 　　-0.317855 　　2.000000 0.400054 -0.320087 -0.275466 -0.284067 　　-0.243598 　　2.250000 0.329940 -0.244935 -0.212786 -0.218538 　　-0.189482 　　2.500000 0.275895 -0.190295 -0.166841 -0.170744 　　-0.149563 　　2.750000 0.233602 -0.150068 -0.132704 -0.135399 　　-0.119703 　　3.000000 0.200020 -0.120024 -0.106973 -0.108868 　　-0.097048 　　3.250000 0.172989 -0.097256 -0.087300 -0.088657 　　-0.079618 　　3.500000 0.150956 -0.079757 -0.072054 -0.073042 　　-0.066030 　　3.750000 0.132790 -0.066124 -0.060087 -0.060818 　　-0.055305 　　4.000000 0.117655 -0.055371 -0.050580 -0.051129 　　-0.046743 　　4.250000 0.104924 -0.046789 -0.042945 -0.043363 　　-0.039833 　　4.500000 0.094123 -0.039866 -0.036750 -0.037072 　　-0.034202 　　4.750000 0.084885 -0.034226 -0.031675 -0.031926 　　-0.029571 　　xn=5.000000 yn=0.076927 我这还有一个程序，从别的论坛上看到的，你看看对你有没有帮助！ function change_step_RK(fun); % 4阶变步长龙格-库塔法解dy=fun(x,y)型方程 % Example: % fun=inline('cos(x)+sin(y)'); % change_step_RK(fun); AbsTol=1e-4; h=0.01; p=4; % p是对应的阶数 x0=0; % x0是起点 xN=1; % xN是终点 y0=0; % 初值 x=x0; y=y0; n=1; p21=2^p-1; while x(end)AbsTol; x1=x2; y1=y2; h=h/2; [x2,y2]=RK_f(fun,x(n),y(n),h/2); end end [xa,ya]=RK_f(fun,h,x(n),y(n)); x(n+1)=xa; y(n+1)=ya; n=n+1; end plot(x,y,'k'); function [xa,ya]=RK_f(fun,h,x,y); % 4阶龙格-库塔法解下一点的值 k1=fun(x,y); k2=fun(x+h/2,y+h*k1/2); k3=fun(x+h/2,y+h*k2/2); k4=fun(x+h,y+h*k3); xa=x+h; ya=y+h*(k1+k2*2+2*k3+k4)/6; 拉格朗日 function y=lagrange(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end SOR迭代法的Matlab程序 function [x]=SOR_iterative(A,b) % 用SOR迭代求解线性方程组,矩阵A是方阵 x0=zeros(1,length(b)); % 赋初值 tol=10^(-2); % 给定误差界 N=1000; % 给定最大迭代次数 [n,n]=size(A); % 确定矩阵A的阶 w=1; % 给定松弛因子 k=1; % 迭代过程 while k<=N x(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1); for i=2:n x(i)=(1-w)*x0(i)+w*(b(i)-A(i,1:i-1)*x(1:i-1)'-A(i,i+1:n)*x0(i+1:n)')/A(i,i); end if max(abs(x-x0))<=tol fid = fopen('SOR_iter_result.txt', 'wt'); fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n'); fprintf(fid,'迭代次数: %d次\n\n',k); fprintf(fid,'x的值\n\n'); fprintf(fid, '%12.8f \n', x); break; end k=k+1; x0=x; end if k==N+1 fid = fopen('SOR_iter_result.txt', 'wt'); fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n'); fprintf(fid,'迭代次数: %d次\n\n',k); fprintf(fid,'超过最大迭代次数，求解失败！'); fclose(fid); end Matlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。 (Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'） Yn+1=Yn+h*f(Xn,Yn) 另外根据微分中值定理，存在0

                    本文档为【龙格库塔法的编程】，请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑，
                    图片更改请在作品中右键图片并更换，文字修改请直接点击文字进行修改，也可以新增和删除文档中的内容。 
 该文档来自用户分享，如有侵权行为请发邮件ishare@vip.sina.com联系网站客服，我们会及时删除。

                    [版权声明] 本站所有资料为用户分享产生，若发现您的权利被侵害，请联系客服邮件isharekefu@iask.cn，我们尽快处理。

                    本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权，请谨慎使用。

                    网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传，仅限个人学习分享使用，禁止用于任何广告和商用目的。
                

下载需要：免费已有0 人下载

立即下载

龙格库塔法的编程

你可能还喜欢