四阶龙格-库塔法解微分方程(C++)
- 格式:docx
- 大小:52.33 KB
- 文档页数:2


首先,龙格库塔法(Runge-Kutta method)是一种常用的数值积分算法,它可以用来数值地求解一阶常微分方程(ODE)。
对于一维ODE,经典的四阶龙格库塔法(RK4)是最常用的方法之一。
下面是一个使用C++实现经典四阶龙格库塔法解一阶微分方程组的示例代码:cpp#include <iostream>#include <vector>// 定义微分方程组dy/dx = f(x, y)std::vector<double> f(double x, const std::vector<double>& y) {std::vector<double> dy(y.size());// 此处为具体的微分方程组表达式,请根据实际需求来定义dy[0] = x * y[0]; // 示例中假设有一个一维微分方程y' = x*yreturn dy;}// 经典四阶龙格库塔法void rk4(double x0, double y0, double h, int numSteps) {double x = x0;double y = y0;for (int i = 0; i < numSteps; ++i) {std::vector<double> k1 = f(x, std::vector<double>{y});std::vector<double> k2 = f(x + h / 2, std::vector<double>{y + h * k1[0] / 2});std::vector<double> k3 = f(x + h / 2, std::vector<double>{y + h * k2[0] / 2});std::vector<double> k4 = f(x + h, std::vector<double>{y + h * k3[0]});y += h * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6;x += h;std::cout << "x: " << x << ", y: " << y << std::endl;}}int main() {double x0 = 0.0; // 初值xdouble y0 = 1.0; // 初值ydouble h = 0.1; // 步长int numSteps = 10; // 迭代次数rk4(x0, y0, h, numSteps);return 0;}请注意,这只是一个简单的示例代码,用于演示如何使用经典四阶龙格库塔法求解一维微分方程组。
数值分析期末试题及答案试题一:1. 简答题(共10分)a) 什么是数值分析?它的主要应用领域是什么?b) 请简要解释迭代法和直接法在数值计算中的区别。
2. 填空题(共10分)a) 欧拉方法是一种______型的数值解法。
b) 二分法是一种______法则。
c) 梯形法则是一种______型的数值积分方法。
3. 计算题(共80分)将以下函数进行数值求解:a) 通过使用二分法求解方程 f(x) = x^3 - 4x - 9 = 0 的近似解。
b) 利用欧拉方法求解微分方程 dy/dx = x^2 + 2x + 1, y(0) = 1 在 x = 1 处的解。
c) 使用梯形法则计算积分∫[0, π/4] sin(x) dx 的近似值。
试题二:1. 简答题(共10分)a) 请解释什么是舍入误差,并描述它在数值计算中的影响。
b) 请解释牛顿插值多项式的概念及其应用。
2. 填空题(共10分)a) 数值稳定性通过______号检查。
b) 龙格-库塔法是一种______计算方法。
c) 零点的迭代法在本质上是将方程______转化为______方程。
3. 计算题(共80分)使用牛顿插值多项式进行以下计算:a) 已知插值节点 (-2, 1), (-1, 1), (0, 2), (1, 4),求在 x = 0.5 处的插值多项式值。
b) 已知插值节点 (0, 1), (1, 2), (3, 7),求插值多项式,并计算在 x = 2 处的值。
c) 使用 4 阶龙格-库塔法求解微分方程 dy/dx = x^2 + 1, y(0) = 1。
答案:试题一:1. a) 数值分析是研究使用数值方法解决数学问题的一门学科。
它的主要应用领域包括数值微积分、数值代数、插值和逼近、求解非线性方程、数值积分和数值解微分方程等。
b) 迭代法和直接法是数值计算中常用的两种方法。
迭代法通过反复迭代逼近解,直到满足所需精度为止;而直接法则通过一系列代数运算直接得到解。
//状态方程/X=AX+BU化为微分方程求解#include <stdio.h>#include <conio.h>#include <math.h>#include <stdlib.h>//定义运算步数;//定义步长;#define N 1000000#define h 0.000001double x[9]={1,2,3,4,5,6,7,8,9};double u[4]={0};//定义微分方程:double fx(double x[],double dx,int i){//矩阵A和BdoubleA[9][9]={{0,0,0,0,0,0,1,0.0023,0.0556},{0,0,0,0,0,0,0,0.9991,-0.0421},{0,0,0,0,0 ,0,0,0.0422,1.0007},{0.0001,-32.1478,0,-0.0197,0.004,0.0142,0.9141,1.7506,0.217} ,{32.1193,-0.0754,0,-0.0044,-0.0323,0.0046,-1.7475,0.9664,0.3464},{-1.3556,-1.78 91,0,0.0141,0.0031,-0.2634,0.1234,0.1345,-2.2448},{0,0,0,-0.0099,-0.031,0.0017,-3.2757,0.721,0.2452},{0,0,0,0.0044,-0.002,-0.0007,-0.1167,-0.6853,-0.0207},{0,0, 0,0.0009,0.0101,-0.0008,0.3469,-0.0784,-0.2875}};doubleB[9][4]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0.0541,-0.0044,0.0218,-0.0007},{-0.0028, -0.0362,0.007,-0.0423},{0.0019,-0.0011,-0.3784,0},{0.001,-0.044,0.0154,-0.0339}, {-0.0136,0.001,-0.0013,-0.0001},{0,-0.004,-0.0162,0.0245}};double a[9] = {0.0};double b[9] = {0.0};for(int j=0;j<9;j++){a[i] += A[i][j]*(x[j]+dx);}for(int k=0;k<4;k++){b[i] += B[k][0]*u[k];}return a[i]+b[i] ;}//主函数void main(){double integer[3]={0};double Kx[9][4];int count=0;int n,S;FILE *fp1,*fp,*fp_integer;fp=fopen("sjy.txt","w");fp1=fopen("sjykxy.txt","w");fp_integer = fopen("integer.txt","w");fprintf(fp1,"x1\tx2\tx3\tx4\tx5\tx6\tx7\tx8\tx9\n");fprintf(fp_integer,"S1\tS2\tS3\n");if(fp==NULL||fp1==NULL){printf("Failed to open file.\n");getch();return;}printf("Input the value of const u1,u2,u3,u4(seperated by ,eg 0.1,0.2,0.3,0.4):"); scanf("%f,%f,%f,%f",&u[0],&u[1],&u[2],&u[3]);printf("Input the value of Steps to get different values of xt,yt(S):");scanf("%d",&S);for(int j=1;j<N;++j){printf("%.3lf",j*h);fprintf(fp,"%.3lf",j*h);//四阶龙格库塔法:for(int i=0; i<9; i++){Kx[i][0]=fx(x,0,i);Kx[i][1]=fx(x,h/2*Kx[i][0],i);Kx[i][2]=fx(x,h/2*Kx[i][1],i);Kx[i][3]=fx(x,h*Kx[i][2],i);x[i]=x[i]+(Kx[i][0]+2*Kx[i][1]+2*Kx[i][2]+Kx[i][3])/6*h;}for(i=0;i<9;++i){//三位有效数字printf("\t%.3lf",x[i]);fprintf(fp,"\t%.3lf",x[i]);}printf("\n");fprintf(fp,"\n");count ++;if(count<=S){integer[0] += x[3]*h;integer[1] += x[4]*h;integer[2] += x[5]*h;for(int m =0;m<3;m++)//六位有效数字fprintf(fp_integer,"\t%.6lf",integer[m]);fprintf(fp_integer,"\n");}//取第S步,即时间为S*h的x1,x2,x3,y1,y2,y3随k值的变化;while(j==S){for(n=0;n<9;++n){//三位有效数字fprintf(fp1,"\t%.3lf",x[n]);}break;}continue;}fclose(fp_integer);fclose(fp1);fclose(fp);printf("\nFinished.\nDatas have been saved to \"sjy.txt,sjykxy.txt,integer.txt\".\n");getch();}建议对飞机追踪器(无人机的黑匣子)进行研究,考虑参考无线电定向技术。
四阶龙格库塔法解微分方程一、四阶龙格库塔法解一阶微分方程①一阶微分方程:cos y t ,初始值y(0)=0,求解区间[0 10]。
MATLAB 程序:%%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程%%%%%%%%%%% y'=cost%%%%%%%%%%% y(0)=0, 0≤t ≤10,h=0.01%%%%%%%%%%% y=sinth=0.01;hf=10;t=0:h:hf;y=zeros(1,length(t));y(1)=0;F=@(t,y)(cos(t));for i=1:(length(t)-1)k1=F(t(i),y(i));k2=F(t(i)+h/2,y(i)+k1*h/2);k3=F(t(i)+h/2,y(i)+k2*h/2);k4=F(t(i)+h,y(i)+k3*h);y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4)*h;endsubplot(211)plot(t,y,'-')xlabel('t');ylabel('y');title('Approximation');span=[0,10];p=y(1);[t1,y1]=ode45(F,span,p);subplot(212)plot(t1,y1)xlabel('t');ylabel('y');title('Exact');图1②一阶微分方程:()22*/x t x x t =- ,初始值x(1)=2,求解区间[1 3]。
MATLAB 程序: %%%%%%%%%%% 四阶龙哥库塔法解微分方程%%%%%%%%%%% x'(t)=(t*x-x^2)/t^2%%%%%%%%%%% x(1)=2, 1≤t ≤3, h=1/128%%%%%%%%%%% 精确解:x(t)=t/(0.5+lnt)h=1/128; %%%%% 步长tf=3;t=1:h:tf;x=zeros(1,length(t));x(1)=2; %%%%% 初始值F_tx=@(t,x)(t.*x-x.^2)./t.^2;for i=1:(length(t)-1)k_1=F_tx(t(i),x(i));k_2=F_tx(t(i)+0.5*h,x(i)+0.5*h*k_1);k_3=F_tx((t(i)+0.5*h),(x(i)+0.5*h*k_2));k_4=F_tx((t(i)+h),(x(i)+k_3*h));x(i+1)=x(i)+(1/6)*(k_1+2*k_2+2*k_3+k_4)*h; endsubplot(211)plot(t,x,'-');xlabel('t');ylabel('x');legend('Approximation');%%%%%%%%%%%%%%%%%%%%%%%%%%%% ode45求精确解t0=t(1);x0=x(1);xspan=[t0 tf];[x_ode45,y_ode45]=ode45(F_tx,xspan,x0);subplot(212)plot(x_ode45,y_ode45,'--');xlabel('t');ylabel('x');legend('Exact');图2二、四阶龙格库塔法解二阶微分方程①二阶微分方程:cos y t ,初始值y(0)=0,y'(0)=-1,求解区间[0 10]。