下载文档原格式
姓名：樊元君学号：02 日期：一、实验目的掌握MATLAB语言、C/C++语言编写计算程序的方法、掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。
掌握使用MATLAB程序求解常微分方程问题的方法。
：二、实验内容1、分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。
实验中以下列数据验证程序的正确性。
求，步长h=。
*2、实验注意事项的精确解为，通过调整步长，观察结果的精度的变化^)三、程序流程图：●改进欧拉格式流程图：~|●四阶龙格库塔流程图：]四、源程序：●改进后欧拉格式程序源代码：function [] = GJOL(h,x0,y0,X,Y)format longh=input('h=');…x0=input('x0=');y0=input('y0=');disp('输入的范围是：');X=input('X=');Y=input('Y=');n=round((Y-X)/h);\i=1;x1=0;yp=0;yc=0;for i=1:1:nx1=x0+h;yp=y0+h*(-x0*(y0)^2);%yp=y0+h*(y0-2*x0/y0);%·yc=y0+h*(-x1*(yp)^2);%yc=y0+h*(yp-2*x1/yp);%y1=(yp+yc)/2;x0=x1;y0=y1;y=2/(1+x0^2);%y=sqrt(1+2*x0);%fprintf('结果=%.3f,%.8f,%.8f\n',x1,y1,y);:endend●四阶龙格库塔程序源代码：function [] = LGKT(h,x0,y0,X,Y)。
format longh=input('h=');x0=input('x0=');y0=input('y0=');disp('输入的范围是：');"X=input('X=');Y=input('Y=');n=round((Y-X)/h);i=1;x1=0;k1=0;k2=0;k3=0;k4=0;for i=1:1:n~x1=x0+h;k1=-x0*y0^2;%k1=y0-2*x0/y0;%k2=(-(x0+h/2)*(y0+h/2*k1)^2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);% k3=(-(x0+h/2)*(y0+h/2*k2)^2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);% k4=(-(x1)*(y0+h*k3)^2);%k4=(y0+h*k3)-2*(x1)/(y0+h*k3);%…y1=y0+h/6*(k1+2*k2+2*k3+k4);%y1=y0+h/6*(k1+2*k2+2*k3+k4);%x0=x1;y0=y1;y=2/(1+x0^2);%y=sqrt(1+2*x0);%fprintf('结果=%.3f,%.7f,%.7f\n',x1,y1,y);end·end*五、运行结果：改进欧拉格式结果：;}四阶龙格库塔结果：步长分别为：和时，不同结果显示验证了步长减少，对于精度的提高起到很大作用，有效数字位数明显增加。
龙格库塔法解微分方程组引言微分方程组是数学中经常遇到的问题，在物理、工程和自然科学中都有广泛的应用。
为了求解微分方程组，我们需要利用数值方法来逼近解析解。
本文将介绍一种常用的数值方法——龙格库塔法（Runge-Kutta method），并探讨如何利用该方法来解微分方程组。
龙格库塔法概述龙格库塔法是一种迭代数值方法，用于求解常微分方程的初值问题。
它的主要思想是将微分方程的解进行离散化，将其转化为一系列的逼近值。
龙格库塔法的基本步骤如下：1.确定步长h和迭代次数n。
2.初始化初始条件，并假设第一个逼近值为y(xi)。
3.依次计算每个逼近值，直到得到y(xi+n*h)为止。
在每个迭代步骤中，龙格库塔法根据前一步的逼近值来计算下一步的逼近值。
该方法具有高阶精度和较好的稳定性，在实际应用中广泛使用。
单一微分方程的龙格库塔法首先，我们来看如何利用龙格库塔法来解一阶常微分方程。
以方程dy/dx = f(x, y)为例，其中f(x, y)为给定的函数。
步骤一：确定步长和迭代次数选择合适的步长h和迭代次数n来进行数值计算。
步长h决定了离散化的精度，而迭代次数n决定了逼近解的数目。
步骤二：初始化条件并计算逼近值设初始条件为y(x0) = y0，其中x0为起始点，y0为起始点处的函数值。
我们先通过欧拉法计算出y(x0 + h)的逼近值，然后再通过该逼近值来计算下一个逼近值。
逼近值的计算公式如下：k1 = h * f(x0, y0)k2 = h * f(x0 + h/2, y0 + k1/2)k3 = h * f(x0 + h/2, y0 + k2/2)k4 = h * f(x0 + h, y0 + k3)y(x0 + h) = y0 + 1/6 * (k1 + 2k2 + 2k3 + k4)步骤三：重复步骤二直到得到y(xi+n*h)依次利用上一步计算出的逼近值来计算下一个逼近值，直到得到y(xi+n*h)为止。
微分方程组的龙格库塔法对于一阶微分方程组的初值问题，我们可以将其转化为向量形式。
首先，龙格库塔法（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;}请注意，这只是一个简单的示例代码，用于演示如何使用经典四阶龙格库塔法求解一维微分方程组。
四阶龙格库塔法(runge-kutta)求振动方程四阶龙格库塔法可以用于求解振动方程。
振动方程通常是一个二阶微分方程，可以转化为一个一阶微分方程组来使用四阶龙格库塔法求解。
假设我们要求解的振动方程为：m * d^2x/dt^2 + c * dx/dt + k * x = 0其中，m是质量，c是阻尼系数，k是弹性系数，x是位移，t 是时间。
首先，将振动方程转化为一阶微分方程组。
引入新的变量v来表示速度，则有：dx/dt = vdv/dt = (-c * v - k * x) / m接下来，使用四阶龙格库塔法求解该一阶微分方程组。
首先将时间区间分割成n个小区间，假设每个小区间的宽度为h。
则利用四阶龙格库塔法的公式可以得到如下更新规则：k1 = h * f(t, x, v)k2 = h * f(t + h/2, x + k1/2, v + l1/2)k3 = h * f(t + h/2, x + k2/2, v + l2/2)k4 = h * f(t + h, x + k3, v + l3)其中，f(t, x, v)表示的是微分方程组的右侧。
根据上述更新规则，可以得到新的位移和速度：x_new = x + (k1 + 2k2 + 2k3 + k4) / 6v_new = v + (l1 + 2l2 + 2l3 + l4) / 6重复以上步骤，可以一步一步地求解出整个时间区间内的位移和速度。
需要注意的是，初始条件x(0)和v(0)需要提供。
总结起来，使用四阶龙格库塔法求解振动方程的步骤如下：1. 将二阶微分方程转化为一阶微分方程组。
2. 设置初始条件x(0)和v(0)。
3. 将时间区间分割成n个小区间，计算每个小区间的宽度h。
4. 根据上述更新规则，逐步求解出整个时间区间内的位移和速度。
希望以上解答对您有所帮助！。
四阶龙格库塔法解微分方程一、四阶龙格库塔法解一阶微分方程①一阶微分方程：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]。
4阶经典龙格库塔公式求解微分方程4阶龙格库塔法（也称为四阶Runge-Kutta方法）是一种常用于求解常微分方程的数值方法。
它是由德国数学家卡尔·龙格以及他的学生马丁·威尔海姆·库塔于1901年独立提出的。
以下是详细的介绍。
1.问题描述我们考虑一个典型的常微分方程初值问题：dy/dx = f(x, y)，y(x0) = y0其中，x0是给定的初始点，y(x)是我们要求解的函数，f(x,y)是已知的函数。
2.方法原理四阶龙格库塔方法的基本思想是通过使用全区间的函数信息来估计下一步的函数值。
在每个步骤中，我们计算四个不同的估计量，然后将它们组合起来以获取最终的解。
具体来说，我们首先计算在初始点x0上的斜率k1=f(x0,y0)。
然后我们计算在x0+h/2处的斜率k2=f(x0+h/2,y0+h/2*k1)，其中h是步长。
以此类推，我们计算k3和k4分别在x0+h/2和x0+h处的斜率。
然后，我们通过加权组合这些斜率来计算下一个点的函数值y1：y1=y0+(h/6)*(k1+2*k2+2*k3+k4)。
3.算法步骤以下是使用四阶龙格库塔法求解常微分方程的详细步骤：步骤1：给定初始条件 y(x0) = y0，选择步长 h 和积分终点 x_end。
步骤2：计算积分步数n：n = (x_end - x0)/h。
步骤3：设置变量x=x0和y=y0，并将步长设置为h。
步骤4：对每个步数i=1,2,...,n，执行以下计算：4.1：计算斜率k1=f(x,y)。
4.2：计算斜率k2=f(x+h/2,y+h/2*k1)。
4.3：计算斜率k3=f(x+h/2,y+h/2*k2)。
4.4：计算斜率k4=f(x+h,y+h*k3)。
4.5：计算下一个点的函数值y1=y+(h/6)*(k1+2*k2+2*k3+k4)。
4.6：将x更新为x+h，y更新为y1步骤5：重复步骤4，直到达到积分终点 x_end。
四阶龙格库塔法（Runge-Kutta ）求解微分方程张晓颖（天津大学 材料学院 学号：1012208027）1 引言计算传热学中通常需要求解常微分方程。
这类问题的简单形式如下：{),(')(00y x f y y x y == （1）虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中的多数微分方程需要采用数值解法求解。
初值问题（1）的数值解法有个基本特点，它们采取“步进式”，即求解过程顺着节点排序一步一步向前推进。
这类算法是要给出用已知信息y n 、 y n −i ……计算y n +1的递推公式即可。
2 龙格库塔法（Runge-Kutta ）介绍假设对于初值问题（1）有解 y = y (x ) ，用 Taylor 展开有：......)(!3)(!2)()()(321+'''+''+'+=+n n n n n x y h x y h x y h x y x y （2） 龙格库塔法（Runge-Kutta ）实质上是间接的使用 Taylor 级数法的一种方法。
对于差商hx y x y n n )()(1-+，根据微分中值定理，存在 0 < θ < 1 ，使得：)()()(1h x y hx y x y n n θ+'=-+ （3）于是对于 y = y (x ) ，可得到：))(,()()(1h x y h x hf x y x y n n n n θθ+++=+ （4）设))(,(*h x y h x f K n n θθ++=，为区间 [x n , x n +1 ] 上的平均斜率。
四阶龙格库塔格式中的*K 由下式计算得到：⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+),()2,2()2,2(),()22(6342312143211hK y h x f K K h y h x f K K h y h x f K y x f K K K K K h y y n n n n nn n n n n （5） 四阶龙格库塔法（Runge-Kutta ）的每一步需要四次计算函数值f ，其截断误差更低，计算的精度更高。
四阶龙格库塔法求解微分方程组题目matlab 在数学和工程学中，微分方程组求解是一个常见的问题。
其中，四阶龙格库塔法是一种常用的求解微分方程组的数值方法。
在本文中，我们将介绍如何使用 MATLAB 实现四阶龙格库塔法求解微分方程组。
首先，我们需要定义微分方程组。
假设我们要求解以下微分方程组：y1' = -2*y1 + y2y2' = -y1 + 2*y2可以将这个微分方程组写成向量形式：Y' = f(Y)其中，Y = [y1; y2]f(Y) = [-2*y1 + y2; -y1 + 2*y2]接下来，我们需要编写四阶龙格库塔法的代码。
具体步骤如下：1. 定义初始值 Y0 和时间间隔 dt。
2. 编写 f(Y) 的函数代码。
3. 使用四阶龙格库塔法计算下一个时间步的 Y 值。
4. 重复步骤 3 直到计算到指定的时间点。
以下是完整的 MATLAB 代码：% 定义初始值 Y0 和时间间隔 dtt0 = 0;tf = 10;dt = 0.01;Y0 = [1; 0];% 编写 f(Y) 的函数代码function dy = f(Y)dy = [-2*Y(1) + Y(2); -Y(1) + 2*Y(2)];end% 使用四阶龙格库塔法计算下一个时间步的 Y 值function Y = rk4(Y, f, dt)k1 = dt*f(Y);k2 = dt*f(Y + 0.5*k1);k3 = dt*f(Y + 0.5*k2);k4 = dt*f(Y + k3);Y = Y + (k1 + 2*k2 + 2*k3 + k4)/6;end% 计算 Y 的值Y = Y0;for t = t0:dt:tfY = rk4(Y, @f, dt);disp(Y);end以上代码将输出从 t0 到 tf 时间段内每个时间步的 Y 值。
您可以根据需要修改微分方程组和时间间隔等参数。
《MATLAB语言及应用》大作业姓名:学号:学院:班级:题目编号:2013年10月134阶Runge-Kutta法求解一阶常微分方程。
一、Runge-Kutta法的数学理论龙格-库塔(Runge-Kutta)方法就是一种在工程上应用广泛的高精度单步算法。
由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。
该算法就是构建在数学支持的基础之上的。
龙格库塔方法的理论基础来源于泰勒公式与使用斜率近似表达微分,它在积分区间多预计算出几个点的斜率,然后进行加权平均,用做下一点的依据,从而构造出了精度更高的数值积分计算方法。
如果预先求两个点的斜率就就是二阶龙格库塔法,如果预先取四个点就就是四阶龙格库塔法。
一阶常微分方程可以写作:y'=f(x,y),使用差分概念。
(Yn+1-Yn)/h= f(Xn,Yn)推出(近似等于,极限为Yn')Yn+1=Yn+h*f(Xn,Yn)另外根据微分中值定理,存在0<t<1,使得Yn+1=Yn+h*f(Xn+th,Y(Xn+th))这里K＝f(Xn+th,Y(Xn+th))称为平均斜率,龙格库塔方法就就是求得K的一种算法。
利用这样的原理,经过复杂的数学推导(过于繁琐省略),可以得出截断误差为O(h^5)的四阶龙格库塔公式:K1＝f(Xn,Yn);K2=f(Xn+h/2,Yn+(h/2)*K1);K3=f(Xn+h/2,Yn+(h/2)*K2);K4=f(Xn+h,Yn+h*K3);Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6)二、Runge-Kutta的算法与流程图在龙格-库塔法中,四阶龙格-库塔法的局部截断误差约为0(h5),被广泛应用于解微分方程的初值问题。
其算法公式为:[]()()()()112341213243/622,/2,/2*,/2,/2*,,n n n n n n n n n n y y h k k k k k f x y k f x h y h k k f x h y h k k f x h y hk +⎧⎫=++++⎪⎪=⎪⎪⎪⎪=++⎨⎬⎪⎪=++⎪⎪⎪⎪=++⎩⎭流程图:(1)、四阶龙格-库塔方法流程图:(2)、实例求解流程图:三、Runge-Kutta的Matlab实现function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次就是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数)n=floor((b-a)/h);%求步数x(1)=a;%时间起点y(:,1)=y0;%赋初值,可以就是向量,但就是要注意维数for ii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii));k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);k4=ufunc(x(ii)+h,y(:,ii)+h*k3);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龙格库塔方法进行数值求解end四、Runge-Kutta的算例实现例:求解常微分方程d(y)/d(x)=-2*y+2*x*x+2*x ,0≤x≤0、5 ,y(0)=1、编辑m文件Untitled3:fun=inline('-2*y+2*x*x+2*x');[x,y]=ode45(fun,[0,0、5],1)>> Untitled3x =0、01250、02500、03750、0625 0、0750 0、0875 0、1000 0、1125 0、1250 0、1375 0、1500 0、1625 0、1750 0、1875 0、2000 0、2125 0、2250 0、2375 0、2500 0、2625 0、2750 0、2875 0、3000 0、31250、33750、35000、36250、37500、38750、40000、41250、42500、43750、45000、46250、47500、48750、5000 y =1、00000、97550、95190、9073 0、8864 0、8663 0、8471 0、8287 0、8112 0、7944 0、7785 0、7633 0、7489 0、7353 0、7224 0、7103 0、6989 0、6883 0、6783 0、6690 0、6605 0、6526 0、6454 0、63880、6277 0、6231 0、6191 0、6157 0、6130 0、6109 0、6093 0、6084 0、6080 0、6083 0、6091 0、6104 0、6124 0、6148 0、6179。
四阶龙格——库塔法2013-2014（1）专业课程实践论文题目：四阶龙格—库塔法2'''()11()()()......()()2!!pp p n n n n n h h y y x hy x y x y x O h p ++=+++++ 进行离散化，所得计算公式必为p 阶方法，式中'''''()(,),()(,)(,)(,)....x y x f x y y x f x y f x y f x y ==++由此，我们能够想到，通过提高泰勒展开式的阶数，可以得到高精度的数值方法，从理论上讲，只要微分方程的解()y x 充分光滑，泰勒展开方法可以构造任意的有限阶的计算公式，但事实上，具体构造这种公式往往相当困难，因为符合函数(,())f x y x 的高阶导数常常是很烦琐的，因此，泰勒展开方法一般不直接使用，但是我们可以间接使用泰勒展开方法，求得高精度的计算方法。
首先，我们对欧拉公式和改进欧拉公式的形式作进一步的分析。
如果将欧拉公式和改进的欧拉公式改写成如下的形式：欧拉公式{111(,)n n n n y y hK K f x y +==+改进的欧拉公式11211()22n n y y h K K +=++， 1(,)n n K f x y =，21(,)n n K f x h y hK =++。
这两组公式都是用函数(,)f x y 在某些点上的值的线性组合来计算1()n y x +的近似值1n y +，欧拉公式每前进一步，就计算一次(,)f x y 的值。
另一方面它是1()n y x +在n x 处的一阶泰勒展开式，因而是一阶方法。
改进的欧拉公式每前进一步，需要计算两次(,)f x y 的值。
另一方面它在(,)n n x y 处的泰勒展开式与1()n y x +在n x 处的泰勒展开式的前三项完全相同，因而是二阶方法。
这启发我们考虑用函数(,)f x y 在若干点上的函数值的线性组合来构造计算公式。
1.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析（1-1），（1-2）（1-3）（1-4）（1-5）（1-6）（1-7）（1-8）（1-9）（1-10）经过循环计算由推得……每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计算。
4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常精准，稳定，且易于编程。
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-2 经典四阶龙格库塔法解一阶微分方程算法程序调试图2.高斯列主元法解线性方程组2.1高斯列主元法解线性方程组算法分析使用伪代码编写高斯消元过程：for k=1 to n-1 dofor i=k+1 to nl<=a(i,k)/a(k,k)for j=k to n doa(i,j)<=a(i,j)-l*a(k,j)end %end of for jb(i)<=b(i)-l*b(k)end %end of for iend %end of for k最后得到A，b可以构成上三角线性方程组接着使用回代法求解上三角线性方程组因为高斯消元要求a（k,k）≠0（k=1,2,3……n-1）这就需要对高斯消元过程进行完善，即使用高斯列主元法：其步骤为：①找主元：计算,并记录其所在行r ,②交换第r行与第k行；③以第k行为工具行处理以下各行，使得从第k列的第k+1行到第n行的元素全部为0；④得到增广矩阵的上三角线性方程组；⑤使用回代法对上三角线性方程组进行求解2.2高斯列主元法解线性方程组流程图图2-1 高斯列主元法解线性方程组流程图2.3高斯列主元法解线性方程组程序代码#include<iostream>#include <cmath>#define N 3using namespace std;void main(){int i,j,k,n,p;float t,s,m,a[N][N],b[N],x[N];cout<<"请输入方程组的系数"<<endl;for(i=0;i<N;i++){for(j=0;j<N;j++)cin>>a[i][j];}cout<<"请输入方程组右端的常数项:"<<endl;for(i=0;i<N;i++)cin>>b[i];for(j=0;j<N-1;j++){ p=j;for(i=j+1;i<N;i++) //寻找系数矩阵每列的最大值{if(fabs(a[i][j])>fabs(a[p][j]))p=i;}if(p!=j) //交换第p行与第j行{for(k=0;k<N;k++){t=a[j][k];a[j][k]=a[p][k];a[p][k]=t;} //交换常数项的第p行与第j行t=b[p];b[p]=b[j];b[j]=t;} //把系数矩阵第j列a[j][j]下面的元素变为0for(i=j+1;i<N;i++){ m=-a[i][j]/a[j][j];for(n=j;n<N;n++)a[i][n]=a[i][n]+a[j][n]*m;b[i]=b[i]+b[j]*m;}} //求方程组的一个解x[N-1]=b[N-1]/a[N-1][N-1]; //回代法求方程组其他解for(i=N-2;i>=0;i--){s=0.0;for(j=N-1;j>i;j--){s=s+a[i][j]*x[j];x[i]=(b[i]-s)/a[i][i];}}cout<<"方程组的解如下:"<<endl;for(i=0;i<=N-1;i++)cout<<x[i]<<endl;}2.4高斯列主元法解线性方程组程序调试结果图示：求解下列方程组图2-2 高斯列主元法解线性方程组程序算法调试图3.牛顿法解非线性方程组3.1牛顿法解非线性方程组算法说明牛顿法主要思想是用进行迭代。
因此首先需要算出的雅可比矩阵，再求过求出它的逆，当它达到某个精度时即停止迭代。
具体算法如下：首先设已知。
①：计算函数（3-1）②：计算雅可比矩阵（3-2）（3-3）③：求线性方程组的解。
④：计算下一点重复上述过程。
3.2牛顿法解非线性方程组流程图图3-1 牛顿法解非线性方程组流程图3.3牛顿法解非线性方程组程序代码#include<iostream>#include<cmath>#define N 2 // 非线性方程组中方程个数、未知量个数#define Epsilon 0.0001 // 差向量1范数的上限#define Max 100 //最大迭代次数using namespace std;const int N2=2*N;int main(){void ff(float xx[N],float yy[N]);//计算向量函数的因变量向量yy[N]void ffjacobian(float xx[N],float yy[N][N]);//计算雅克比矩阵yy[N][N] void inv_jacobian(float yy[N][N],float inv[N][N]);//计算雅克比矩阵的逆矩阵invvoid newdundiedai(float x0[N], float inv[N][N],float y0[N],float x1[N]);//由近似解向量 x0 计算近似解向量 x1floatx0[N]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errorno rm;int i,j,iter=0;//如果取消对x0的初始化，撤销下面两行的注释符，就可以由键盘向x0读入初始近似解向量//for( i=0;i<N;i++)// cin>>x0[i];cout<<"初始近似解向量："<<endl;for (i=0;i<N;i++)cout<<x0[i]<<" ";cout<<endl;cout<<endl;do{iter=iter+1;cout<<"第 "<<iter<<" 次迭代开始"<<endl;//计算向量函数的因变量向量 y0ff(x0,y0);//计算雅克比矩阵 jacobianffjacobian(x0,jacobian);//计算雅克比矩阵的逆矩阵 invjacobianinv_jacobian(jacobian,invjacobian);//由近似解向量 x0 计算近似解向量 x1newdundiedai(x0, invjacobian,y0,x1);//计算差向量的1范数errornormerrornorm=0;for (i=0;i<N;i++)errornorm=errornorm+fabs(x1[i]-x0[i]);if (errornorm<Epsilon) break;for (i=0;i<N;i++)x0[i]=x1[i];} while (iter<Max);return 0;}void ff(float xx[N],float yy[N]){float x,y;int i;x=xx[0];y=xx[1];yy[0]=x*x-2*x-y+0.5;yy[1]=x*x+4*y*y-4;cout<<"向量函数的因变量向量是： "<<endl;for( i=0;i<N;i++)cout<<yy[i]<<" ";cout<<endl;cout<<endl;}void ffjacobian(float xx[N],float yy[N][N]){float x,y;int i,j;x=xx[0];y=xx[1];//jacobian have n*n elementyy[0][0]=2*x-2;yy[0][1]=-1;yy[1][0]=2*x;yy[1][1]=8*y;cout<<"雅克比矩阵是： "<<endl;for( i=0;i<N;i++){for(j=0;j<N;j++)cout<<yy[i][j]<<" ";cout<<endl;}cout<<endl;}void inv_jacobian(float yy[N][N],float inv[N][N]) {float aug[N][N2],L;int i,j,k;cout<<"开始计算雅克比矩阵的逆矩阵："<<endl;for (i=0;i<N;i++){ for(j=0;j<N;j++)aug[i][j]=yy[i][j];for(j=N;j<N2;j++)if(j==i+N) aug[i][j]=1;else aug[i][j]=0;}for (i=0;i<N;i++){ for(j=0;j<N2;j++)cout<<aug[i][j]<<" ";cout<<endl;}cout<<endl;for (i=0;i<N;i++){for (k=i+1;k<N;k++){L=-aug[k][i]/aug[i][i];for(j=i;j<N2;j++)aug[k][j]=aug[k][j]+L*aug[i][j];}}for (i=0;i<N;i++){ for(j=0;j<N2;j++)cout<<aug[i][j]<<" ";cout<<endl;}cout<<endl;for (i=N-1;i>0;i--){for (k=i-1;k>=0;k--){L=-aug[k][i]/aug[i][i];for(j=N2-1;j>=0;j--)aug[k][j]=aug[k][j]+L*aug[i][j];}}for (i=0;i<N;i++){ for(j=0;j<N2;j++)cout<<aug[i][j]<<" ";cout<<endl;}cout<<endl;for (i=N-1;i>=0;i--)for(j=N2-1;j>=0;j--)aug[i][j]=aug[i][j]/aug[i][i];for (i=0;i<N;i++){ for(j=0;j<N2;j++)cout<<aug[i][j]<<" ";cout<<endl;for(j=N;j<N2;j++)inv[i][j-N]=aug[i][j];}cout<<endl;cout<<"雅克比矩阵的逆矩阵： "<<endl;for (i=0;i<N;i++){ for(j=0;j<N;j++)cout<<inv[i][j]<<" ";cout<<endl;}cout<<endl;}void newdundiedai(float x0[N], float inv[N][N],float y0[N],float x1[N]) {int i,j;float sum=0;for(i=0;i<N;i++){ sum=0;for(j=0;j<N;j++)sum=sum+inv[i][j]*y0[j];x1[i]=x0[i]-sum;}cout<<"近似解向量："<<endl;for (i=0;i<N;i++)cout<<x1[i]<<" ";cout<<endl;cout<<endl;}3.4牛顿法解非线性方程组程序调试结果图示：图3-2 牛顿法解非线性方程组程序算法调试图图3-3 牛顿法解非线性方程组程序算法调试图图3-4 牛顿法解非线性方程组程序算法调试图4.龙贝格求积分算法4.1龙贝格求积分算法分析生成j>=k 的近似积分结果逼近表，并以R(j+1,j+1)为最终解来逼近积分。
