经典四阶龙格库塔法解一阶微分方程组
MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程
姓名：樊元君学号：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)为止。
微分方程组的龙格库塔法对于一阶微分方程组的初值问题，我们可以将其转化为向量形式。
c++经典四阶龙格库塔法解一阶微分方程组
首先，龙格库塔法（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)求振动方程
四阶龙格库塔法(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阶经典龙格库塔公式求解微分方程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)求解微分方程
四阶龙格库塔法（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 在数学和工程学中，微分方程组求解是一个常见的问题。
其中，四阶龙格库塔法是一种常用的求解微分方程组的数值方法。
在本文中，我们将介绍如何使用 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 值。
您可以根据需要修改微分方程组和时间间隔等参数。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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-2 经典四阶龙格库塔法解一阶微分方程算法程序调试图⎪⎪⎩⎪⎪⎨⎧+=+=yx dtdy y x dt dx232⎩⎨⎧==4)0(6)0(y x2.高斯列主元法解线性方程组2.1高斯列主元法解线性方程组算法分析使用伪代码编写高斯消元过程：for k=1 to n-1 do for 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 j b(i)<=b(i)-l*b(k) end %end of for i end %end of for k最后得到A ，b 可以构成上三角线性方程组 接着使用回代法求解上三角线性方程组因为高斯消元要求a （k,k ）≠0（k=1,2,3……n-1）这就需要对高斯消元过程进行完善，即使用高斯列主元法：其步骤为：①找主元：计算 (){}1max ||()k pk a p k n -≤≤,并记录其所在行r ,||rk a = (){}1max ||()k pk a p k n -≤≤②交换第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]下面的元素变为0 for(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高斯列主元法解线性方程组程序调试结果图示：求解下列方程组341-11182522321321321=++=-+=++x x x x x x x x x图2-2 高斯列主元法解线性方程组程序算法调试图3.牛顿法解非线性方程组3.1牛顿法解非线性方程组算法说明牛顿法主要思想是用(1)()()1()()()k k k k x x F x F x +-'=- (0,1,...).k = 进行迭代 。
因此首先需要算出()F x 的雅可比矩阵()F x '，再求过()F x '求出它的逆1()F x -'，当它达到某个精度时即停止迭代。
具体算法如下：首先设k P 已知。
①：计算函数12(,)()(,)k k k k k f p q F P f p q ⎡⎤=⎢⎥⎣⎦②：计算雅可比矩阵1122(,)(,)()(,)(,)k k k k k k k k k f p q f p q x yJ P f p q f p q x y ∂∂⎡⎤⎢⎥∂∂⎢⎥=∂∂⎢⎥⎢⎥∂∂⎣⎦③：求线性方程组 ()()k k J P P F P ∆=-的解P ∆。
④：计算下一点1k k P P P +=+∆重复上述过程。
（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],errornorm;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)为最终解来逼近积分。