龙格库塔法 C语言程序

合集下载

四阶龙格库塔求解边界层问题(C语言)

四阶龙格库塔求解边界层问题(C语言)

题目: 设52/ 1.110/m s μρ-=⨯的气体以10/v m s ∞=的速度以零攻角定常扰流长度为L =1m 的大平板,用数值解讨论边界层内的流动规律。

如图所示,沿板面方向取x 坐标,板的法线方向取y 坐标。

在流体力学中,介绍了存在相似性解的二维平面不可压缩流体的边界层方程:C.E :0=∂∂+∂∂yx u υM.E : 221y udx dp y u x u u ∂∂+-=∂∂+∂∂νρυ0p y∂=∂ 边界条件为:y = 0;u = v = 0y =∞;u =v ∞ ( u = u (x) 为x 点处壁面势流流速 )在外边界上221122e e p v p v c ρρ∞+=+=所以 00e dpdx+=对于平板扰流,则平板扰流的边界层方程可简化为 C.E :0=∂∂+∂∂yx u υ M.E : 22u u uu x y yυν∂∂∂+=∂∂∂ 其定解的边界条件为y = 0;u = v = 0y =∞;u =v ∞为了求解边界层方程,我们可以利用“相似性解”的方法,对其进行转化,由C.E ,引进流函数),(y x ψ,令xy u ∂∂-=∂∂=ψυψ, 由M.E 式的偏导阶次,可望减少自变量数目令()yg x η==()()f f v g x η∞== 其中2x x ηη∂=-∂1()y g x η∂=∂ 由()v g x f ψ∞=,得()'v g x fu v f y yψ∞∞∂∂===∂∂ ()()(')2v g x f v g x v f f x x xψη∞∞∂∂=-==-∂∂ 所以,(')"2v u v f f x x x η∞∞∂∂==-∂∂(')"()v u v f f y y g x ∞∞∂∂==∂∂222(")"'()()v v u f f y y g x g x ∞∞∂∂==∂∂ 将其都代入M.E ,化简可得"'"0f ff += 这就使原方程变化为:"'"0f ff +=此时的边界条件为:η = 0;f(0) = 0 , f’(0) = 0η = ∞;f’(∞) = 1那么,如何利用编辑程序的方法求解这个变化后的边界层微分方程呢?一. 解方程的基本思路为了简化运算,此时边界层微分方程化简成:"'"0f ff +=,边界条件不变。

龙格库塔算法

龙格库塔算法

龙格库塔算法龙格库塔算法(Runge-Kutta method)是一种常用的数值解微分方程的方法,其基本原理是通过逐步逼近的方式,根据初始条件和微分方程的表达式,计算出方程的近似解。

该方法具有较高的精度和稳定性,在科学计算、物理模拟、工程建模等领域得到广泛应用。

龙格库塔算法的核心思想是将微分方程的解按照一定的步长进行离散化,从而将连续的求解问题转化为离散的迭代过程。

具体来说,龙格库塔算法通过计算函数在一定步长内的平均斜率,来估计下一个点的函数值。

这个平均斜率是通过多次计算函数在不同点上的导数得到的,从而提高了计算的精度。

龙格库塔算法的一般形式可以表示为:k1 = f(tn, yn)k2 = f(tn + h/2, yn + h/2 * k1)k3 = f(tn + h/2, yn + h/2 * k2)k4 = f(tn + h, yn + h * k3)yn+1 = yn + h/6 * (k1 + 2k2 + 2k3 + k4)其中,tn是当前时间点,yn是当前函数值,h是步长,f是微分方程的表达式。

通过多次迭代,可以逐渐逼近微分方程的解。

龙格库塔算法的优点在于其精确度较高,可以通过调整步长来控制计算的精度和效率。

此外,该算法具有较好的数值稳定性,可以有效处理非线性、刚性或高阶微分方程等复杂问题。

因此,在科学和工程计算中,龙格库塔算法被广泛应用于各种数值模拟和求解问题。

需要注意的是,龙格库塔算法并非万能的,对于一些特殊的问题,可能存在数值不稳定性或计算精度不够的情况。

此外,算法的步长选择也需要根据具体问题进行调整,过小的步长会增加计算量,而过大的步长可能导致精度下降。

因此,在使用龙格库塔算法时,需要根据具体问题的特点和要求来选择合适的步长和算法参数,以获得满意的计算结果。

总结起来,龙格库塔算法是一种常用的数值解微分方程的方法,具有较高的精度和稳定性。

通过离散化和迭代的方式,可以逐步逼近微分方程的解。

龙格-库塔法

龙格-库塔法

四阶龙格-库塔法求解常微分方程的初值问题1.算法原理对于一阶常微分方程组的初值问题⎪⎪⎪⎩⎪⎪⎪⎨⎧=⋯⋯==⋯⋯=⋯⋯⋯⋯=⋯⋯=0020********'212'2211'1)(,,)(,)())(,),(),(,()())(,),(),(,()())(,),(),(,()(n n n n n n n y x y y x y y x y x y x y x y x f x y x y x y x y x f x y x y x y x y x f x y , 其中b x a ≤≤。

若记Tn Tn Tn y x f y x f y x f y x f y y y y x y x y x y y x y )),(,),,(),,((),(),,,())(),(),(()(2102010021⋯⋯=⋯⋯=⋯⋯=,,则可将微分方程组写成向量形式⎩⎨⎧=≤≤=0')()),(,()(y a y b x a x y x f x y微分方程组初值问题在形式上和单个微分方程处置问题完全相同,只是数量函数在此变成了向量函数。

因此建立的单个一阶微分方程初值问题的数值解法,可以完全平移到求解一阶微分方程组的初值问题中,只不过是将单个方程中的函数转向向量函数即可。

标准4阶R-K 法的向量形式如下:⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧++=++=++==++++=+),()21,2()21,2(),()22(61342312143211K y h x hf K K y h x hf K K y h x hf K y x hf K K K K K y y n n n n n n n n n n 其分量形式为n j K y K y K y h x hf K K y K y K y h x hf K K y K y K y h x hf K y y y x hf K K K K K y y n ni i i i j j n nii i i j j n nii i i j j ni i i i j j j j j j i j i j ,,2,1).,,,;(),2,2,2;2(),2,2,2;2(),,,,;(),22(6132321314222212131212111221143211,1,⋯⋯=⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧+⋯⋯+++=+⋯⋯+++=+⋯⋯+++=⋯⋯=++++=++,,2.程序框图3.源代码%该函数为四阶龙格-库塔法function [x,y]=method(df,xspan,y0,h)%df为常微分方程,xspan为取值区间,y0为初值向量,h为步长x=xspan(1):h:xspan(2);m=length(y0);n=length(x);y=zeros(m,n);y(:,1)=y0(:);for i=1:n-1k1=feval(df,x(i),y(:,i));k2=feval(df,x(i)+h/2,y(:,i)+h*k1/2);k3=feval(df,x(i)+h/2,y(:,i)+h*k2/2);k4=feval(df,x(i)+h,y(:,i)+h*k3);y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;end%习题9.2clear;xspan=[0,1];%取值区间h=0.05;%步长y0=[-1,3,2];%初值df=@(x,y)[y(2);y(3);y(3)+y(2)-y(1)+2*x-3];[xt,y]=method(df,xspan,y0,h)syms t;yp=t*exp(t)+2*t-1;%微分方程的解析解yp1=xt.*exp(xt)+2*xt-1%计算区间内取值点上的精确解[xt',y(1,:)',yp1']%y(1,:)为数值解,yp1为精确解ezplot(yp,[0,1]);%画出解析解的图像hold on;plot(xt,y(1,:),'r');%画出数值解的图像4.计算结果。

四阶龙格_库塔算法的C语言实现_毋玉芝

四阶龙格_库塔算法的C语言实现_毋玉芝

四阶龙格—库塔算法的C 语言实现毋玉芝(焦作财会学校)摘 要 本文叙述了四阶龙格—库塔算法的C 语言实现过程、数据存储及其结果的曲线显示,并以具体实例说明了这一过程。

关键词 龙格 库塔算法 数据存储 曲线显示在科学技术中常常需要求解常微分方程的定解问题,这就需要一种合适的数值解法求出常微分方程的解。

在诸多数值算法中龙格—库塔算法具有较高的精确度,是一种优先选取的算法。

1.四阶龙格—库塔算法简述龙格—库塔方法实际上是间接地使用台劳级数法的一种技术。

龙格—库塔算法的数学描述如下: y n +1=y n +h ·(K 1+2K 2+2K 3+K 4)/6; K 1=f (x n ,y n ); K 2=f (x n +h /2,y n +K 1·h /2); K 3=f (x n +h /2,y n +K 2·h /2); K 4=f (x n +h ,y n +K 3·h );其中: h 表示计算过程中选取的步长;K 1表示x n 点处的斜率; K 2表示利用K 1求得的(x n +h /2)点处的斜率; K 3表示利用K 2求得的(x n +h /2)点处的斜率; K 4表示利用K 3求得的(x n +h )点处的斜率;2.C 语言的实现过程2.1 算法实现龙格—库塔算法关键是选择步长h ,必须根据题目的要求选出合适的步长,这对龙格—库塔算法结果的精确度及其平滑性尤为重要。

在选择了恰当的步长后,利用上述迭代表达式,并根据题目要求的迭代次数,或求解的精度,利用C 语言加以实现。

2.2 数据存储由于此计算结果数据庞大,程序运行后数据不可能一屏显示,鉴于此首先利用fopen ()函数创建并打开一文本文件,利用fprintf ()函数将数据存储到数据文件,可将此文件打印输出,以检验结果的正确性。

实现过程如下:if (fp ==N ULL ){printf (″Can ′t open this file n ″);exit (0);}for (i =0;i <=j ;i ++){t =ts +i *h ;2001年3月第1期 焦作大学学报JO U RNA L OF JIAOZ UO U NI VERSIT Y №.1Mar .200156 焦 作 大 学 学 报 2001年3月fprintf(fp,″t=%3.2f n″,t);fprintf(fp,″x=%8.5f ″,x[i]);fprintf(fp,″y=%8.5f ″,y[i]);fprintf(fp,″z=%8.5f n″,z[i]);}fclose(fp);2.3 曲线显示根据上述所取得的结果,并根据坐标的适当变换即可以描绘出结果的平滑曲线。

三阶Runge-Kutta方法

三阶Runge-Kutta方法

2012-2013(1)专业课程实践论文三阶Runge-Kutta方法王曹旭,0818180106,R数学08-1班常微分方程初值问题的数值解法是求方程(1)的解在点列1(0,1,)n n n x x h n -=+= 上的近似值n y ,这里n h 是1n x -到n x 的步长,一般略去下标记为h 。

00(,)()dyf x y dxy x y⎧=⎪⎨⎪=⎩(1)三阶Runge-Kutta 方法,它的计算公式是:⎪⎪⎩⎪⎪⎨⎧+++=++==+++=+))((,(),,(),,(),(2131213322111sK rK qh y qh x f K phK y ph x f K y x f K K K K h y y n n n n n n n n λλλ R K-方法的优点是:单步法、精度高,计算过程便于改变步长,缺点是计算量较大,每前进一步需要计算三次函数值f 。

在用龙格库塔方法时,要注意N 的选择要合适,N 太大,会使计算量加大,N 太小,h 较大,可能会使误差增大。

因此选择合适的N 很重要。

我们要在考虑精度的基础上,选择合适的N 。

在此,用c 语言实现了三阶Runge-Kutta 方法。

三阶Runge-Kutta流程图#include<stdlib.h>#include<stdio.h>/*n表示几等分,n+1表示他输出的个数*/int RungeKutta(double y0,double a,double b,int n,double *x,double *y,double (*function)(double,double)){double h=(b-a)/n,k1,k2,k3;int i;// x=(double*)malloc((n+1)*sizeof(double));// y=(double*)malloc((n+1)*sizeof(double));x[0]=a;y[0]=y0;for(i=0;i<n;i++){x[i+1]=x[i]+h;k1=function(x[i],y[i]);k2=function(x[i]+h/2,y[i]+h*k1/2);k3=function(x[i]+h,y[i]-h*k1+2*h*k2);y[i+1]=y[i]+h*(k1+4*k2+k3)/6;}return 1;}double function(double x,double y){return y-2*x/y;}//例子求y'=y-2*x/y(0<x<1);y0=1;int main(){ int i;double x[7],y[7];printf("用三阶龙格-库塔方法\n");RungeKutta(1,0,1,6,x,y,function);for(i=0;i<7;i++)printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);return 1;}例1.求解方程⎩⎨⎧=<<+='1,10,0y x y x y 步长2.0=h解:(1)将程序序中rerurn y-2*x/y 改为rerurn x+y(2)输入,,,,0n b a y 的值:;5;1;0;1(3)显示输出结果:435019.3y 650070.2043616.2583310.1242667.1000000.1000000.1800000.0600000.0400000.0200000.0000000.05432054320============y y y y y x x x x x x 时,时,时,时,时,时,例2.求解方程⎩⎨⎧=<<+='110),1/(30y x x y y 步长2.0=h解:(1)将程序序中rerurn y-2*x/y 改为rerurn 3*y/(1+x)(2) 输入n b a y ,,,0的值:;5;1;0;1(3) 显示输出结果:961755.7805800.5079391.4734787.2724242.1000000.1000000.1800000.0600000.0x 400000.0x 200000.0000000.0543210543210============y y y y y y x x x x 时,时,时,时,时,时,。

龙格-库塔方法(Runge-Kutta)

龙格-库塔方法(Runge-Kutta)

龙格-库塔方法〔Runge-Kutta〕3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的一般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值方法,若用显式单步法(3.2.2)当,即数值求积用左矩形公式,它就是Euler法(3.1.2),方法只有一阶精度,若取(3.2.3)就是改进Euler法,这时数值求积公式是梯形公式的一种近似,计算时要用二个右端函数f的值,但方法是二阶精度的.若要得到更高阶的公式,则求积分时必须用更多的f值,根据数值积分公式,可将(3.2.1)右端积分表示为注意,右端f中还不能直接得到,需要像改进Euler法(3.1.11)一样,用前面已算得的f值表示为(3.2.3),一般情况可将(3.2.2)的 表示为(3.2.4)其中这里均为待定常数,公式(3.2.2),(3.2.4)称为r级的显式Runge-Kutta法,简称R-K方法.它每步计算r个f值(即由前面(i-1)个已算出的表示,故公式是显式的.例),而ki如当r=2时,公式可表示为(3.2.5)其中.改进Euler 法(3.1.11)就是一个二级显式R-K 方法.参数取不同的值,可得到不同公式.3.2.2 二、三级显式R-K 方法对r=2的显式R-K 方法(3.2.5),要求选择参数,使公式的精度阶p 尽量高,由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6)令,对(3.2.6)式在处按Taylor 公式展开,由于将上述结果代入(3.2.6)得要使公式(3.2.5)具有的阶p=2,即,必须(3.2.7)即由此三式求的解不唯一.因r=2,由〔3.2.5〕式可知,于是有解(3.2.8)它表明使(3.2.5)具有二阶的方法很多,只要都可得到二阶精度R-K方法.若取,则,则得改进Euler法(3.1.11),若取,则得,此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)与中点公式(3.2.9)是两个常用的二级R-K方法,注意二级R-K方法只能达到二阶,而不可能达到三阶.因为r=2只有4个参数,要达到p=3则在(3.2.6)的展开式中要增加3项,即增加三个方程,加上(3.2.7)的三个方程,共计六个方程求4个待定参数,验证得出是无解的.当然r=2,p=2的R-K方法(3.2.5)当取其他数时,也可得到其他公式,但系数较复杂,一般不再给出.对r=3的情形,要计算三个k值,即其中将按二元函数在处按Taylor公式展开,然后代入局部截断误差表达式,可得可得三阶方法,其系数共有8个,所应满足的方程为(3.2.10)这是8个未知数6个方程的方程组,解也是不唯一的,通常.一种常见的三级三阶R-K方法是下面的三级Kutta方法:(3.2.11)附:R-K 的三级Kutta 方法程序如下function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满足c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K 方法与步长的自动选择利用二元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K 方法,其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =,2221(,(,))n n n n k f x a h y b hf x y =++,而33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定,Taylor 展开分析局部截断误差,使得精度达到四阶,即误差为5()O h 。

c求解四阶龙格库塔法

c求解四阶龙格库塔法

某一地区的病菌传染,三种人员的人数的状态方程,即可能受传染的人数x 1,已被传染的病的人数x2,及已治愈的人数x3,并假设他们是因接触而传染。

设:α是x1中单位时间内的传染系数,β是x2中单位时间内被治愈的比例系数,则得以下状态方程:f`(x1)=-α*x1*x2;f`(x2)=α*x1*x2-β*x2;f`(x3)=β*x2;其中f`(x)是加速度本题中,α=0.01,β=0.027初始条件:t=0, x1=620人, x2=10人, x3=70人选步长:h=0.025, 用四阶RK法模拟。

-------------程序清单-------------------#include <stdio.h>#include <conio.h>#define A 0.01#define B 0.027void myFun(double *x,double *f){f[0]=-A*x[0]*x[1];f[1]=A*x[0]*x[1]-B*x[1];f[2]=B*x[1] ;}int rht4(double *rs,double h,double x[],int n,int Num){int j,i,ii;double ps[4]={0,0.5,0.5,1};double bs[4]={1.0/6,1.0/3,1.0/3,1.0/6};double K[3];double y[3];double f[3];for(i=0;i<n;i++){y[i]=x[i];rs[i]=y[i];}for(j=1;j<Num;j++){for(ii=0;ii<4;ii++){ /*龙格-库塔*/for(i=0;i<n;i++){x[i]=rs[(j-1)*n+i]+ps[ii]*f[i]*h; /*特别注意rs[(j-1)*n+i] 而不是x[i] 不是前次的*/}myFun(x,f);/*算每次的加速度*/for(i=0;i<n;i++){y[i]+=f[i]*h*bs[ii];}} /*end 龙格-库塔*/for(i=0;i<n;i++){x[i]=y[i];rs[j*n+i]=y[i];}}}int main(){double t=0,x[3]={620.0,10.0,70.0},h=0.025,*rh;int Num,i,j;double f[3];printf("Input Num: ");scanf("%d",&Num);rh=(double *)malloc(Num*3*sizeof(double));rht4(rh,h,x,3,Num) ;printf("x1 x2 x3-----------------------------------------------------"); for(i=0;i<Num;i++){printf(" ");for(j=0;j<3;j++)printf("%f ",rh[i*3+j]);}getch();return 0;}----------end---------------。

龙格库塔算法解微分方程组c语言

龙格库塔算法解微分方程组c语言

/***************************************************************************This program is to solve the initial value problem of following systemof differential equations:dx/dt=x+2*y,x(0)=0,dy/dt=2*x+y,y(0)=2,x(0.2) and y(0。

2) are to be calculated****************************************************************************/#include<stdio。

h>#include<math。

h〉#define steplength 0。

1 //步?长¡è可¨¦根¨´据Y需¨¨要°a调Ì¡Â整?;#define FuncNumber 2 //FuncNumber为a微¡é分¤?方¤?程¨¬的Ì?数ºy目?;void main(){float x[200],Yn[20][200],reachpoint;int i;x[0]=0;Yn[0][0]=0;Yn[1][0]=2; //初?值¦Ì条¬?件t;reachpoint=0。

2; //所¨´求¨®点Ì?可¨¦根¨´据Y需¨¨要°a调Ì¡Â整?;void rightfunctions(float x ,float *Auxiliary,float *Func);void Runge_Kutta(float *x,float reachpoint, float(*Yn)[200]);Runge_Kutta(x ,reachpoint, Yn);printf("x ");for(i=0;i<=(reachpoint—x[0])/steplength;i++)printf("%f ",x[i]);printf("\nY1 ”);for(i=0;i〈=(reachpoint—x[0])/steplength;i++)printf("%f ”,Yn[0][i]);printf(”\nY2 ”);for(i=0;i〈=(reachpoint—x[0])/steplength;i++)printf("%f ”,Yn[1][i]);getchar();}void rightfunctions(float x ,float *Auxiliary,float *Func)//当Ì¡À右®¨°方¤?程¨¬改?变À?时º¡À,ê?需¨¨要°a改?变À?;{Func[0]=Auxiliary[0]+2*Auxiliary[1];Func[1]=2*Auxiliary[0]+Auxiliary[1];}void Runge_Kutta(float *x,float reachpoint, float(*Yn)[200]){int i,j;float Func[FuncNumber],K[FuncNumber][4],Auxiliary[FuncNumber];for(i=0;i<=(reachpoint—x[0])/steplength;i++){for(j=0;j<FuncNumber;j++)Auxiliary[j]=*(Yn[j]+i);rightfunctions(x[i],Auxiliary,Func);for(j=0;j<FuncNumber;j++){K[j][0]=Func[j];Auxiliary[j]=*(Yn[j]+i)+0.5*steplength*K[j][0];}rightfunctions(x[i],Auxiliary,Func);for(j=0;j<FuncNumber;j++){K[j][1]=Func[j];Auxiliary[j]=*(Yn[j]+i)+0.5*steplength*K[j][1];}rightfunctions(x[i],Auxiliary,Func);for(j=0;j〈FuncNumber;j++){K[j][2]=Func[j];Auxiliary[j]=*(Yn[j]+i)+steplength*K[j][2];}rightfunctions(x[i],Auxiliary,Func);for(j=0;j〈FuncNumber;j++)K[j][3]=Func[j];for(j=0;j〈FuncNumber;j++)Yn[j][i+1]=Yn[j][i]+(K[j][0]+2*K[j][1]+2*K[j][2]+K[j][3])*steplength/6.0;x[i+1]=x[i]+steplength;}}。

龙格库塔法的c

龙格库塔法的c

龙格库塔法的c++编程CODE:#include<stdlib.h>#include<stdio.h>/*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<n;i++){x[i+1]=x[i]+h;k1=function(x[i],y[i]);k2=function(x[i]+h/2,y[i]+h*k1/2);y[i+1]=y[i]+h*k2;}break;case 3:for(i=0;i<n;i++){x[i+1]=x[i]+h;k1=function(x[i],y[i]);k2=function(x[i]+h/2,y[i]+h*k1/2);k3=function(x[i]+h,y[i]-h*k1+2*h*k2);y[i+1]=y[i]+h*(k1+4*k2+k3)/6;}break;case 4:for(i=0;i<n;i++){x[i+1]=x[i]+h;k1=function(x[i],y[i]);k2=function(x[i]+h/2,y[i]+h*k1/2);k3=function(x[i]+h/2,y[i]+h*k2/2);k4=function(x[i]+h,y[i]+h*k3);y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;}break;default:return 0;}return 1;}double function(double x,double y){return y-2*x/y;}//例子求y'=y-2*x/y(0<x<1);y0=1;/*int main(){double x[6],y[6];printf("用二阶龙格-库塔方法\n");RungeKutta(1,0,1,5,x,y,2,function);for(int i=0;i<6;i++)printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);printf("用三阶龙格-库塔方法\n");RungeKutta(1,0,1,5,x,y,3,function);for(i=0;i<6;i++)printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);printf("用四阶龙格-库塔方法\n");RungeKutta(1,0,1,5,x,y,4,function);for(i=0;i<6;i++)printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);return 1;龙格库塔求解微分方程数值解工程中很多的地方用到龙格库塔求解微分方程的数值解,龙格库塔是很重要的一种方法,尤其是四阶的,精确度相当的高。

使用C语言解常微分方程 C ODE【范本模板】

使用C语言解常微分方程 C ODE【范本模板】

解常微分方程姓名:Vincent年级:2010,学号:1033****,组号:5(小组),4(大组)1. 数值方法:我们的实验目标是解常微分方程,其中包括几类问题。

一阶常微分初值问题,高阶常微分初值问题,常微分方程组初值问题,二阶常微分方程边值问题,二阶线性常微分方程边值问题。

对待上面的几类问题,我们分别使用不同的方法。

• 初值问题使用 龙格-库塔 来处理 • 边值问题用打靶法来处理 • 线性边值问题有限差分法初值问题我们分别使用• 二阶 龙格—库塔 方法 • 4阶 龙格—库塔 方法 来处理一阶常微分方程。

理论如下:对于这样一个方程'()(,)y t f t y =当h 很小时,我们用泰勒展开,12111111(,)(,)(,)k k k k i i k k j j i ij k hf t y k hf t h y k k hf a b a b t h y h k -===++=++∑当我们选择正确的参数 a[i][j],b [i ][j]之后,就可以近似认为这就是泰勒展开. 对于二阶,我们有:()01()()y t h y t h Af Bf+=++其中010(,)(,)P f f t y f f t y h Qhf ==++经过前人的计算机经验,我们有,选择 A=1/2,B=1/2,则 P=1,Q=1,于是又如下形式,这也使休恩方法的表达式。

所以我们称其为 龙格(库塔)休恩方法()()()(,)(,(,))2hy t h y t f t y f t h y hf t y +=++++对于4阶龙格方法,我们有类似的想法, 我们使用前人经验的出的系数,有如下公式112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y k k k k k f t y h h k f t y k h h k f t y k k f t h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩对于高阶微分方程及微分方程组 我们用• 4阶龙格-库塔方法来解对于一个如下的微分方程组'111'1(,,,),(,,,).n nn n y f t y y y f t y y ⎧=⎪⎨⎪=⎩ 101,00,0(),()n n y t y y t y=⎧⎪⎨⎪=⎩ 我们可以认为是一个一阶向量微分方程,所以可以用龙格-库塔方法的向量形式解。

龙格-库塔方法

龙格-库塔方法

8.2 龙格-库塔方法8.2.1 二阶龙格-库塔方法常微分方程初值问题:做在点的泰勒展开:这里。

取,就有(8.11) 截断可得到近似值的计算公式,即欧拉公式:若取,式(8.11)可写成:或(8.12)截断可得到近似值的计算公式:或上式为二阶方法,一般优于一阶的欧拉公式(8.2),但是在计算时,需要计算在点的值,因此,此法不可取。

龙格-库塔设想用在点和值的线性组合逼近式(8.12)的主体,即用(8.13)逼近得到数值公式:(8.14)或更一般地写成对式(8.13)在点泰勒展开得到:将上式与式(8.12)比较,知当满足时有最好的逼近效果,此时式(8.13)-式(8.14)。

这是4个未知数的3个方程,显然方程组有无数组解。

若取,则有二阶龙格-库塔公式,也称为改进欧拉公式:(8.15)若取,则得另一种形式的二阶龙格-库塔公式,也称中点公式:(8.16)从公式建立过程中可看到,二阶龙格-库塔公式的局部截断误差仍为,是二阶精度的计算公式。

类似地,可建立高阶的龙格-库塔公式,同时可知四阶龙格-库塔公式的局部截断误差为,是四阶精度的计算公式。

欧拉法是低精度的方法,适合于方程的解或其导数有间断的情况以及精度要求不高的情况,当解需要高精度时,必须用高阶的龙格-库塔等方法。

四阶龙格-库塔方法应用面较广,具有自动起步和便于改变步长的优点,但计算量比一般方法略大。

为了保证方法的收敛性,有时需要步长取得较小,因此,不适于解病态方程。

8.2.2 四阶龙格-库塔公式下面列出常用的三阶、四阶龙格-库塔计算公式。

三阶龙格-库塔公式(1)(8.17)(2)(8.18)(3)(8.19)四阶龙格-库塔公式(1)(8.20)(2)(8.21)例8.3用四阶龙格-库塔公式(8.20)解初值问题:解:取步长,计算公式为:计算结果列表8.3中。

表8.3 计算结果8.2.3 步长的自适应欧拉方法和龙格-库塔方法在计算时仅用到前一步的值,我们称这样的方法为单步法。

龙格-库塔(Runge-Kutta)方法

龙格-库塔(Runge-Kutta)方法
证明:
( )
dy dy y′ = = f, y′′ = f x + f y = f x + ffy = F dx dx F F dy F F y′′′ = + = +f x y dx x y F = f xx + f x f y + ffyx = f xx + f x f y + ffxy x F 2 2 f = f(fxy + f y f y + ffyy ) = ffxy + ffy + f f yy y
y ( x) = e

x2

x
0
e dt
t2

x
0
e dt 难以求积
t2
ODE数值解的基本思想和方法特点 数值解的基本思想和方法特点
1. 离散化 级数、 用Taylor级数、数值积分和差商逼近导数, 级数 数值积分和差商逼近导数, 将 ODE转化为离散的代数方程 称差分方程 。 转化为离散的代数方程(称差分方程 转化为离散的代数方程 称差分方程)。
(ha2 ) + (ha3 ) +
f 2 2! x
2 2 2
(ha f ) +
2
2
K3 f + ha3 f x + h (a3 b32 ) f + b32 K2 f y 2! 2! + h2a3 (a3 b32 ) f + b32 K2 f xy f xx + h (a3 b32 ) f + b32 K2
2
Euler法 后退 法 ym+1 = ym + hK2 + O(h2 ) K2 = f ( xm + h, ym+1 )

经典四阶龙格库塔法解一阶微分方程组

经典四阶龙格库塔法解一阶微分方程组

1.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析),,(1k k k y x t f f =, )2,2,2(112g hy f h x h t f f k k k +++=)2,2,2(223g hy f h x h t f f k k k +++=),,(334hg y hf x h t f f k k k +++=),,(1k k k y x t g g =)2,2,2(112g hy f h x h t g g k k k +++=)2,2,2(223g hy f h x h t g g k k k +++=),,(334hg y hf x h t g g k k k +++=)22(6)22(64321143211g g g g hy y f f f f hx x k k k k ++++=++++=++1k k t t h +=+经过循环计算由 推得 ……每个龙格-库塔方法都是由一个合适的泰勒方法推导而来,使得其最终全局误差为()N O h ,一种折中方法是每次进行若干次函数求值,从而省去高阶导数计算。

4阶龙格-库塔方法(RK4)是最常用的,它适用于一般的应用,因为它非常精准,稳定,且易于编程。

000,,t x y ()()111222,,,,t x y t x y (1-1)(1-2)(1-3)(1-4)(1-5)(1-6)(1-7)(1-8)(1-9)(1-10)1.2经典四阶龙格库塔法解一阶微分方程流程图图1-1 经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微分方程程序代码:#include <iostream>#include <iomanip>using namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h) {double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;t0=initial[0];x0=initial[1];y0=initial[2];f1=f(t0,x0,y0); g1=g(t0,x0,y0);f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2); g2=g(t0+h/2, x0+h*f1/2,y0+h*g1/2);f3=f(t0+h/2, x0+h*f2/2,y0+h*g2/2); g3=g(t0+h/2,x0+h*f2/2,y0+h*g2/2);f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3); x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6; resu[0]=t0+h;resu[1]=x1;resu[2]=y1;}int main(){double f(double t,double x, double y);double g(double t,double x, double y);double initial[3],resu[3];double a,b,H;double t,step;int i;cout<<"输入所求微分方程组的初值t0,x0,y0:";cin>>initial[0]>>initial[1]>>initial[2];cout<<"输入所求微分方程组的微分区间[a,b]:";cin>>a>>b;cout<<"输入所求微分方程组所分解子区间的个数step:";cin>>step;cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision( 10);H=(b-a)/step;cout<< initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl; for(i=0;i<step;i++){ RK4( f,g ,initial, resu,H);cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<endl;initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];}return(0);}double f(double t,double x, double y){double dx;dx=x+2*y; return(dx); }double g(double t,double x, double y) {double dy; dy=3*x+2*y; return(dy); }1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示:应用所编写程序计算所给例题:其中初值为求解区间为[0,0.2]。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
y1=y1+h*((K1/3)+(2*K2/3));
x1=a+i*h;
x[i]=x1;
printf("x[%d]=%f\t",i,x[i]);
fprintf(fp,"x[%d]=%f\t",i,x[i]);
y[i]=y1;
printf("y[%d]=%f\t",i,y[i]);
fscanf(fp,"%d",n);
printf("\n计算步数:");
printf("n=%d\n",*n) ;
fclose(fp);
}
void RK(int n,float y0,float a,float b,float x[],float y[])
printf("y(xn)\t\t");
fprintf(fp,"y(xn)\t\t");
printf("|y(xn)-yn|\n\n");
fprintf(fp,"|y(xn)-yn|\n\n");
printf("x[0]=%f\t",x[0]) ;
fprintf(fp,"x[0]=%f\t",x[0]) ;
if((fp=fopen(filename,"w"))==NULL)
{
printf("can not open file.\n");
getchar();
exit(0);
}
h=0.1;
x1=a;
y1=y0;
x[0]=x1;
y[0]=y1;
yy[0]=y0;
dy[0]=0;
void RK(int n,float y0,float a,float b,float x[],float y[]);
float f(float x,float y);
float af(float x);
void main()
{
int n;
float a,b,y0,*x,*y;
{
int i;
float x1,y1,h;
float K1,K2;
float *yy,*dy;
FILE *fp;
char *filename="runge_kutta_output.txt";
yy=(float*)malloc((n+1)*sizeof(float));
dy=(float*)malloc((n+1)*sizeof(float));
{
printf("cannot open file runge_kutta.txt!\n");
}
printf("\n二阶Runge_Kutta法求初值:\n");
printf("\n区间:");
fscanf(fp,"%f",a);
fscanf(fp,"%f",b);
printf("[a,b]=[%f,%f]",*a,*b) ;
printf("\n") ;
fscanf(fp,"%f",y0);
printf("\n初值:");
printf("y0=%f\n",*y0) ;
fprintf(fp,"dy(x[0])=%f\n",dy[0]) ;
for(i=1;i<=n;i++)
{
/*K1=f(x1,y1);
K2=f((x1+h),(y1+h*K1));
y1=y1+h*((K1/2)+(K2/2)); */
K1=f(x1,y1);
K2=f((x1+3*h/4),(y1+3*h*K1/4));
fprintf(fp,"y[%d]=%f\t",i,y[i]);
yy[i]=af(x[i]);
printf("y(x[%d])=%f\t",i,yy[i]);
fprintf(fp,"y(x[%d])=%f\t",i,yy[i]);
dy[i]=fabs(y[i]-yy[i]) ;
fun=y-2*x/y;
return fun;
}
float af(float x)
{
float fun;
fun=sqrt(1+2*x);
return fun;
}
输入文件runge_kutta_input.txt中内容:
0
1
1
10
以下为代码部分:
#include<stdio.h>
#include "stdlib.h"
#include<math.h>
#include <conio.h>
void size(int *n) ;
void input(float *a,float *b,float *y0,int *n);
size(&n) ;
x=(float*)malloc((n+1)*sizeof(float));
y=(float*)malloc((n+1)*sizeof(float));
input(&a,&b,&y0,&n);
RK(n,y0,a,b,x,y);
system("pause");
free(x);free(y);
exit(0);
}
fscanf(fp,"%d",n);
fclose(fp);
}
void input(float *a,float *b,float *y0,int *n)
{ FILE *fp;
if((fp=fopen("runge_kutta_input.txt","r"))==NULL)
printf("\n常微分方程y-2*x/y=y'的数值解为:\n");
fprintf(fp,"常微分方程y-2*x/y=y'的数值解为:\n");
printf("\n\txn\t\t");
fprintf(fp,"\n\txn\t\t");
printf("yn\t\t");
fprintf(fp,"yn\t\t");
printf("dy(x[%d])=%f\n",i,dy[i]);
fprintf(fp,"dy(x[%d])=%f\n",i,dy[i]);
}
free(yy);free(dy) ;
fclose(fp);
}
float f(float x,float y)
{
float fun;
printf("y[0]=%f\t",y[0]) ;
fprintf(fp,"y[0]=%f\t",y[0]) ;
printf("y(x[0])=%f\t",yy[0]) ;
fprintf(fp,"y(x[0])=%f\t",yy[0]) ;
printf("dy(x[0])=%f\n",dy[0]) ;
}
void size(int *n) /* 读取维数n */
{
FILE *fp;
if((fp=fopen("runge_kutta_input.txt","r"))==NULL)
{
printf("can not open file.\n");
相关文档
最新文档