精华内容
下载资源
问答
  • 含有隐函数的离散常微分方程求解

    千次阅读 2016-11-28 23:55:24
    微分方程求解问题在大学高等数学中学过,当我碰到这个问题后,我还专门去翻了一下高数的书(都还给老师了ORZ)。微分方程的定义是:凡含有参数,未知函数和未知函数导数 (或微分) 的方程,称为微分方程。而常微分...

      前段时间代码中碰到一个含有隐函数的离散常微分方程的求解问题,在这里和大家讨论一下其求解方法。
      首先,介绍一下什么是常微分方程。微分方程求解问题在大学高等数学中学过,当我碰到这个问题后,我还专门去翻了一下高数的书(都还给老师了ORZ)。微分方程的定义是:凡含有参数,未知函数和未知函数导数 (或微分) 的方程,称为微分方程。而常微分方程是指未知函数是一元函数的微分方程。
      下面是一个简单的常微分方程:
      这里写图片描述
      上式中,f(x)和g(x)都是x的显式函数,a是常数,解上边的微分方程可以得到f(x)的函数表达式。对于连续的情况,即:f(x)和g(x)都是连续函数,我们在大一大二学的微积分或者高数书上都有详细的求解方法。例如:可采用常数变易法,通过常数变易法,可求出一阶线性微分方程的通解。带入初值即可求出其特解。
      常徽分方程是研究自然科学社会科学中的事物运动和变化规律的最基本的数学理论及方法自然界中很多事物的运动规律,如牛顿运动定律,万有引力定律等, 社会科学中很多现象的演变规律如人口发展规律市场均衡价格的变化等 都可用常徽分方程来描述。但是我们在工程中遇到的常微分方程往往不是我们学过的理想形式,下边就分析我遇到的一种情况:
      g(x)是采集的离散值,不知道函数形式,是一个隐函数,总结一下就是:g(x)是隐函数的离散形式的常微分方程求解。对于这个问题,我认为可以用前向欧拉、函数拟合或者构造优化求解模型,利用优化方法进行求解。我这里介绍利用欧拉法和函数拟合的方法求解这个常微分方程。

    1.利用欧拉法求解含有隐函数的微分方程

      这个方法思路非常简单,就是利用差分方程代替微分方程。得到f(x)的每个采样点的值。
      欧拉
      通过上式我们可以得到f(x)在每个点的值。在程序数据处理中,我们得到f(x)的所有值就足够了。下边我用matlab对此方法进行了仿真:

    %% ********用欧拉的方式求解ode *******%%
    %%
    clc;
    clear all;
    close all;
    %% 定义构造函数参量
    func1 = @(x)sin(2*pi*x)+cos(2*pi*x);   
    func2 = @(x)3*x+1;
    func3 = @(x)x.*x+1;
    f1_init = func1(0);           %设置初始点
    f2_init = func2(0);
    f3_init = func3(0);
    constant_A = 2;               %微分方程常系数
    sampling_step = 0.01;         %采样间隔
    % sampling_step = 1;          %采样间隔
    % samp_max = 100;
    samp_max = 10;
    samp_min = 0;
    samp_num = (samp_max-samp_min)/sampling_step;
    sampling_x = linspace(samp_min,samp_max,samp_num);
    %% 通过构造函数获得恒解空间
    % func1_diff = diff(func1,'x');
    % func2_diff = diff(func2,'x');
    % func3_diff = diff(func3,'x');
    func1_diff = @(x)2*pi*cos(2*pi*x)-2*pi*sin(2*pi*x);
    func2_diff = @(x)3;
    func3_diff = @(x)2*x;
    for i=1:samp_num
        g1_fun(i) = func1(sampling_x(i))+constant_A*func1_diff(sampling_x(i));
        g2_fun(i) = func2(sampling_x(i))+constant_A*func2_diff(sampling_x(i));
        g3_fun(i) = func3(sampling_x(i))+constant_A*func3_diff(sampling_x(i));
    end
    % 通过构造的反解
    func1_Inverse = [f1_init];
    func2_Inverse = [f2_init];
    func3_Inverse = [f3_init];
    for k=1:samp_num-1
        func1_Inverse(k+1) = func1_Inverse(k)+(g1_fun(k)-func1_Inverse(k))*sampling_step/constant_A;
        func2_Inverse(k+1) = func2_Inverse(k)+(g2_fun(k)-func2_Inverse(k))*sampling_step/constant_A;
        func3_Inverse(k+1) = func3_Inverse(k)+(g3_fun(k)-func3_Inverse(k))*sampling_step/constant_A;
    end
    %% 画图比较
    raw_fun1_data = func1(sampling_x);   %准确数字
    raw_fun2_data = func2(sampling_x);   %线性可以完全拟合
    raw_fun3_data = func3(sampling_x);
    figure;
    plot(sampling_x,raw_fun1_data,'r');
    hold on;
    plot(sampling_x,func1_Inverse,'b');
    figure;
    plot(sampling_x,raw_fun2_data,'r.');
    hold on;
    plot(sampling_x,func2_Inverse,'b.');
    figure;
    plot(sampling_x,raw_fun3_data,'r.');
    hold on;
    plot(sampling_x,func3_Inverse,'b.');
    

      我测试了正余弦、线性和二次函数,都能很好地完成工作。下边是正余弦函数的测试结果:
    这里写图片描述

    2.利用函数拟合方式求解含有隐函数的微分方程

      由于函数g(x)是离散的隐函数,所以很容易想到可以用函数拟合的方法去得到g(x)的函数形式。我采用了多项式拟合的方法去拟合g(x)。方便地是,matlab提供了函数拟合和求解ode(微分方程)的函数。下边给出代码:

    %% ********用函数拟合的方式求解ode(polyfit()和ode()函数*******%%
    %%
    clc;
    clear all;
    close all;
    %% 定义构造函数参量
    func1 = @(x)sin(2*pi*x)+cos(2*pi*x);   
    func2 = @(x)3*x+1;
    func3 = @(x)x.*x+1;
    f1_init = func1(0);           %设置初始点
    f2_init = func2(0);
    f3_init = func3(0);
    constant_A = 2;               %微分方程常系数
    sampling_step = 0.01;         %采样间隔
    % sampling_step = 1;          %采样间隔
    % samp_max = 100;
    samp_max = 10;
    samp_min = 0;
    samp_num = (samp_max-samp_min)/sampling_step;
    sampling_x = linspace(samp_min,samp_max,samp_num);
    %% 通过构造函数获得恒解空间
    % func1_diff = diff(func1,'x');
    % func2_diff = diff(func2,'x');
    % func3_diff = diff(func3,'x');
    func1_diff = @(x)2*pi*cos(2*pi*x)-2*pi*sin(2*pi*x);
    func2_diff = @(x)3;
    func3_diff = @(x)2*x;
    for i=1:samp_num
        g1_fun(i) = func1(sampling_x(i))+constant_A*func1_diff(sampling_x(i));
        g2_fun(i) = func2(sampling_x(i))+constant_A*func2_diff(sampling_x(i));
        g3_fun(i) = func3(sampling_x(i))+constant_A*func3_diff(sampling_x(i));
    end
    %% 利用多项式拟合g(x)
    g1_fit = polyfit(sampling_x,g1_fun,5);
    g2_fit = polyfit(sampling_x,g2_fun,3);
    g3_fit = polyfit(sampling_x,g3_fun,3);
    g1_func_fit = @(x)g1_fit(1)*x.*x.*x.*x.*x+g1_fit(2)*x.*x.*x.*x+g1_fit(3)*x.*x.*x+g1_fit(4)*x.*x+g1_fit(5)*x+g1_fit(6);
    g2_func_fit = @(x)g2_fit(1)*x.*x.*x+g2_fit(2)*x.*x+g2_fit(3)*x+g2_fit(4);
    g3_func_fit = @(x)g3_fit(1)*x.*x.*x+g3_fit(2)*x.*x+g3_fit(3)*x+g3_fit(4);
    %% ode求解
    % odefun1 = @(x,f1)(g1_fit(1)*x.*x.*x+g1_fit(2)*x.*x+g1_fit(3)*x+g1_fit(4)-f1)/constant_A;
    odefun1 = @(x,f1)(g1_fit(1)*x.*x.*x.*x.*x+g1_fit(2)*x.*x.*x.*x+g1_fit(3)*x.*x.*x+g1_fit(4)*x.*x+g1_fit(5)*x+g1_fit(6))/constant_A;
    odefun2 = @(x,f2)(g2_fit(1)*x.*x.*x+g2_fit(2)*x.*x+g2_fit(3)*x+g2_fit(4)-f2)/constant_A;
    odefun3 = @(x,f3)(g3_fit(1)*x.*x.*x+g3_fit(2)*x.*x+g3_fit(3)*x+g3_fit(4)-f3)/constant_A;
    f1_init = func1(0.01);           %设置初始点
    f2_init = func2(0.01);
    f3_init = func3(0.01);
    [x,f1] = ode45(odefun1,[0.01:0.01: 10],f1_init);
    [x,f2] = ode45(odefun2,[0.01:0.01: 10],f2_init);
    [x,f3] = ode45(odefun3,[0.01:0.01: 10],f3_init);
    %% 画图比较
    sampling_x = [0.01:0.01: 10];
    raw_fun1_data = func1(sampling_x);      %准确数字
    raw_fun2_data = func2(sampling_x);      %线性可以完全拟合
    raw_fun3_data = func3(sampling_x);
    figure;
    plot(sampling_x,raw_fun1_data,'r');
    hold on;
    plot(x,f1,'b');
    figure;
    plot(sampling_x,raw_fun2_data,'r.');
    hold on;
    plot(x,f2,'b.');
    figure;
    plot(sampling_x,raw_fun3_data,'r.');
    hold on;
    plot(x,f3,'b.');
    

      我测试了正余弦、线性和二次函数,都能很好地完成工作。结果会发现,线性,二次函数都能很好的用此方法解决,但是无法完成正余弦函数的求解。因为在函数拟合部分,很难用多项式来拟合正余弦函数。
      具体的实验结果,大家可以用代码跑一下。
      

    展开全文
  • 常微分方程数值求解1. 常微分方程数值求解的一般概念2. 常微分方程数值求解函数3. 常微分方程数值求解函数统一命名格式4. 刚性问题 1. 常微分方程数值求解的一般概念 2. 常微分方程数值求解函数 [t,y]=solver...

    1. 常微分方程数值求解的一般概念

    2. 常微分方程数值求解函数

    [t,y]=solver(filename,tspan,y0,option)
    t:时间向量; y:数值解;
    filename:定义f(t,y)的函数名,该函数必须返回一个列向量;
    tspan形式为[t0,tf],表示求解区间; y0:初始状态向量;
    option:可选参数,用于设置求解属性,常用的属性包括相对误差值RelTol(默认值为10-3和绝对误差值AbsTol(默认值为10-6);

    3. 常微分方程数值求解函数统一命名格式

    odennxx

    • ode:Ordinary Differential Equation,常微分方程
    • nn:数字,代表所用方法的阶数。
      ode23采用2阶龙格-库塔(Runge-Kutta)算法,用3阶公式做误差估计来调节步长,具有低等精度。
      ode45采用4阶龙格-库塔算法,用5阶公式做误差估计来调节步长,具有中等精度。
    • xx:字母,用于标注方法的专门特征。
      ode15s、ode23s中的“s”代表(Stiff),表示函数适用于刚性方程。
      在这里插入图片描述
      e.g.求解微分方程初值问题,并与精确解y(t)=sqrt(t+1)+1,进行比较。
      在这里插入图片描述
    >> f=@(t,y)(y^2-t-2)/4/(t+1);
    >> [t,y]=ode23(f,[0,10],2);
    >> y1=sqrt(t+1)+1;
    >> plot(t,y,'b:',t,y1,'r')
    

    在这里插入图片描述
    e.g.已知一个二阶线性系统的微分方程为:
    在这里插入图片描述

    >> f=@(t,x) [-2,0;0,1]*[x(2);x(1)];
    >> [t,x]=ode45(f,[0,20],[1,0]);
    >> subplot(2,2,1);
    >> plot(t,x(:,2));
    >> subplot(2,2,2);
    >> plot(x(:,2),x(:,1));
    

    在这里插入图片描述

    4. 刚性问题

    有一类常微分方程,其解的分量有的变化很快,有的变化很慢,且相差悬殊,这就是所谓的刚性问题(Stiff)。对于刚性问题,数值解算法必须取很小步长才能获得满意的结果,导致计算量会大大增加。解决刚性问题需要有专门方法。非刚性算法可以求解刚性问题,只不过需要很长的计算时间。
    在这里插入图片描述

    >> lambda=0.01;
    >> f=@(t,y) y^2-y^3;
    >> tic;
    >> [t,y]=ode45(f,[0,2/lambda],lambda);
    >> toc
    时间已过 32.701709 秒。
    >> disp(['ode45计算的点数' num2str(length(t))]);
    ode45计算的点数157
    %tic和toc函数用来记录微分方程求解命令执行的时间,使用tic函数启动计时器,使用toc函数显示从计时器启动到当前所经历的时间。最后还输出计算的点数,运行结果表明这时常微分方程不算很刚性。
    
    >> lambda=1e-5;
    >> f=@(t,y) y^2-y^3;
    >> tic;
    >> f=@(t,y) y^2-y^3;
    >> toc
    时间已过 7.096559 秒。
    >> disp(['ode45计算的点数' num2str(length(t))]);
    ode45计算的点数157
    这时计算时间明显加长,计算的点数剧增,表明这时常微分方程表现为刚性
    
    >> lambda=1e-5;
    >> f=@(t,y) y^2-y^3;
    >> tic;
    >> [t,y]=ode15s(f,[0,2/lambda],lambda);
    >> toc
    时间已过 8.699999 秒。
    >> disp(['ode15s计算的点数' num2str(length(t))]);
    ode15s计算的点数98
    对于刚性问题,选择以“s”结尾的函数,专门用于求解刚性方程。计算时间大大缩短,计算的点数大大减少。
    相同程序,每次的计算时间都不相同
    
    展开全文
  • Matlab 语言编程基础(矩阵的创建及使用方法、矩阵运算基础、线性方程求解、选择、循环{求和、迭代}、函数的定义及使用方法) (40分) 曲线的绘制(掌握一元函数图形的绘制方法,以及参数曲线的绘制方法; 掌握...

    《数学实验》知识点整理

    一、需要掌握的内容

    Matlab 语言编程基础(矩阵的创建及使用方法、矩阵运算基础、线性方程求解、选择、循环{求和、迭代}、函数的定义及使用方法) (40分)

    曲线的绘制(掌握一元函数图形的绘制方法,以及参数曲线的绘制方法;  掌握plot, ezplot,plot3,ezplot3等函数的用法;掌握基本的图形标注命令)  (15分左右)

    曲面的绘制(掌握二元函数图形的绘制方法,以及参数曲面的绘制方法;掌握 mesh,ezmesh等函数的用法)(15分左右)

    非线性方程(组)求解 (掌握用fzero求非线性方程的根,以及用 fsolve求解非线性方程组的方法)(10分)

    数值积分和符号积分(掌握trapz,quadl,dblquad, integral2, integral3等函数的用法) (10分左右)

    常微分方程(组)求解  (掌握利用ode45求解常微分方程和常微分方程组的方法)  (10分)

    二、知识整理

    1. 数值积分数值微分

    1)        常用的数值积分方法

    • 矩形公式
    • 梯形公式
    • Simpson公式
    • Newton-Cotes积分法
    • 高斯求积公式
    • 龙贝格(Romberg)积分法
    • 蒙特卡洛方法(随机模拟的方法)

    2)       Matlab中数值微分和积分的函数:diff

    3)       数值导数

    方向导数gradient

    » clear;x=[1 1.1 1.2 1.3];y=x.^3;

    » dy=diff(y)./diff(x)

    » dy=gradient(y,x)

    » 3*x.^2

    利用梯形法求积分

    • z=trapz(x,y)  其中x 表示积分区间的离散化向量; y   x同维数的向量,表示被积函数;该函数返回积分的近似值
    • 求积分
    • » clear; x=-1:0.1:1; y=exp(-x.^2);
    • » trapz(x,y)

     

    1. 2.        数值积分和符号积分
    • z=trapz(x,y)  其中x 表示积分区间的离散化向量; y是与x同维数的向量,表示被积函数;该函数返回积分的近似值 。
    • 例: 求积分
    • 解 » clear; x=-1:0.1:1; y=exp(-x.^2);
    • » trapz(x,y)
    • >> z=quadl(@(x)exp(-x.^2),-1,1)(int函数)
    • >> z= integral (@(x)exp(-x.^2),-1,1) (可以计算反常积分)
    • 另一个例子:
    • >> fun = @(x) exp(-x.^2).*log(x).^2;
    • >> q = quadl(fun,0,Inf)    % quadl不能求反常积分
    • >> q = integral(fun,0,Inf) % integral能求反常积分
    1. 3.        重积分
    • z=dblquad(Fun,a,b,c,d)  求二元函数 Fun(x,y) 在矩形区域的重积分。
    • z=triplequad(Fun,a,b,c,d,e,f)  求三元函数Fun(x,y,z) 在长方体区域上的三重积分。
    • z=quad2d(Fun,a,b,cx,dx)  求二元函数Fun(x,y)在一般区域上的重积分。a, b为变量x的下、上限;cx, dx为变量y的下、上限函数(自变量为x)。
    • z= integral2 (Fun,a,b,cx,dx) 类似quad2d
    • z= integral3 (Fun,a,b,cx,dx,exy,fxy) 求三元函数Fun(x,y,z)在一般区域上的三重积分 。

    例题:

     

    1. 4.        微分方程

    微分方程:含有未知函数及其某些阶导数以及自变量本身的方程称为微分方程

    常微分方程:未知函数是一元函数

    偏微分方程:含有偏导数的微分方程,其解为多元函数u(t,x,y,z)。

    微分方程组:联系一些未知函数x(t), y(t), z(t),  … 的一组微分方程。

    微分方程的阶:微分方程中出现的未知函数的导数的最高阶数

    n  (2)常系数线性微分方程的特征根法

    n  线性方程:  y(n)  + a1 (t) y(n-1)  + …

    n      + an-1 (t) y’ + an (t) y = b(t)

    n  常系数方程: 若ai (t) (i =1, …,n) 与t无关。

    n  齐次方程: 若b(t)=0。

    n        y(n)  + a1 y(n-1)  + … + an-1 y’ + an y = 0

    n  线性常系数齐次微分方程的解可用特征根法求得.

    n         ln+a1 ln-1+ … +an-1 l+an=0

    n  非齐次方程的解为一个特解和相应的齐次方程通解的叠加。

    n  变系数方程可尝试常数变易法。

    n  例6.1  求x’’+ 0.2 x’+3.92x = 0的通解

    n  解  特征方程为l2 + 0.2l +3.92=0

    n    » roots([1 0.2 3.92]

    n   求得共轭复根 a ±bi=-0.1±1.9774i,

    ode45函数

     

     

    例6.2  解微分方程                                                     y’ -y+2t/y=0, y(0)=1(初值向量1, 0<t<4 (自变量的初值04

    将方程整理为标准形式y’ = y-2t/y

     

    程序:

    odefun= @(t,y)y-2*t/y ;

    >> [t,y]=ode45(odefun,[0,4],1); [t,y]

    >> plot(t,y,'o')

    >> ode45(odefun,0:1:4,1);%(输出结点列向量)

    >> [t,y]=ode45(odefun,0:1:4,1);[t,y]

    例二

     

     

    【0 30】自定义x的范围,【1;0.5】初值解的值

     

     

    例三(高阶)

     

    Y0是初值解

     

    1. 5.        边值问题解法 

     

     

    6.5求解边值问题

    首先改写为标准形式。

         y(1)= z, y(2)= z’, 则方程为

    y’(1)=y(2), y’(2) = -y(2)sin(y(1))

    边界条件为

    ya(1)=0, yb(1)+2=0

    程序:eg6_5.m

     %根据z初始值预估:z=-1,z’=0

    clear;close;

    sinit=bvpinit(0:4,[-1;0])

    %[-1;0]是常数猜测值z=-1, z’=0

    • odefun=@(t,y)[y(2);-y(2)*sin(y(1))];

    bcfun=@(ya,yb)[ya(1);yb(1)+2];

    sol=bvp5c(odefun,bcfun,sinit) %注意sol的域名

    t=linspace(0,4,101);

    y=deval(sol,t);

    plot(t,y(1,:),sol.x,sol.y(1,:),'o',sinit.x,sinit.y(1,:),'s')

    legend('解曲线','解点','粗略解')

     

     

     

     

     

     

    转载于:https://www.cnblogs.com/star-lit-sky-shines/p/10191604.html

    展开全文
  • 这里对使用python求解常微分方程提供两种思路,一种是自己编程实现欧拉法,改进欧拉法或者四阶龙格库塔,这样有助于理解上述三种数值计算方法的原理;一种是调用python已有的库,不再重复造轮子。本文对上述两种思路...

    这里对使用python求解常微分方程提供两种思路,一种是自己编程实现欧拉法,改进欧拉法或者四阶龙格库塔,这样有助于理解上述三种数值计算方法的原理;一种是调用python已有的库,不再重复造轮子。

    本文对上述两种思路都给出代码示例,并进行比较;同时针对单个微分方程和含有多个微分方程的微分方程组给出代码示例

    1. 常微分方程定义

    凡含有参数,未知函数和未知函数导数 (或微分) 的方程,称为微分方程。

    未知函数是一元函数的微分方程称作常微分方程

    未知数是多元函数的微分方程称作偏微分方程。

    微分方程中出现的未知函数最高阶导数的阶数,称为微分方程的阶数。

    2. 调用现有的库

    scipy中提供了用于解常微分方程的函数odeint(),完整的调用形式如下:

    scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, \

    rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0,hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)

    实际使用中,还是主要使用前三个参数,即微分方程的描写函数、初值和需要求解函数值对应的的时间点。接收数组形式。这个函数,要求微分方程必须化为标准形式,即

    ,实际操作中会发现,高阶方程的标准化工作,其实是解微分方程最主要的工作。

    示例1:单个微分方程

    import math

    import numpy as np

    from scipy.integrate import odeint

    import matplotlib.pyplot as plt

    def func(y, t):

    return t * math.sqrt(y)

    YS=odeint(func,y0=1,t=np.arange(0,10.1,0.1))

    t=np.arange(0,10.1,0.1)

    plt.plot(t, YS, label='odeint')

    plt.legend()

    plt.show()

    结果如下:

    图.1

    示例2:微分方程组

    与单个微分方程不同的是,此时的函数变成了向量函数

    from scipy.integrate import odeint

    import numpy as np

    def lorenz(w, t, p, r, b):

    # 给出位置矢量w,和三个参数p, r, b计算出

    # dx/dt, dy/dt, dz/dt的值

    x, y, z = w

    # 直接与lorenz的计算公式对应

    return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])

    t = np.arange(0, 30, 0.01) # 创建时间点

    # 调用ode对lorenz进行求解, 用两个不同的初始值

    track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))

    track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))

    # 绘图

    from mpl_toolkits.mplot3d import Axes3D

    import matplotlib.pyplot as plt

    fig = plt.figure()

    ax = Axes3D(fig)

    ax.plot(track1[:,0], track1[:,1], track1[:,2])

    ax.plot(track2[:,0], track2[:,1], track2[:,2])

    plt.show()

    结果如下

    图.2

    3. 自己编程实现欧拉法/改进欧拉法/四阶龙格库塔

    示例 3:单个函数使用四阶龙格库塔

    import math

    import numpy as np

    import matplotlib.pyplot as plt

    def runge_kutta(y, x, dx, f):

    """ y is the initial value for y

    x is the initial value for x

    dx is the time step in x

    f is derivative of function y(t)

    """

    k1 = dx * f(y, x)

    k2 = dx * f(y + 0.5 * k1, x + 0.5 * dx)

    k3 = dx * f(y + 0.5 * k2, x + 0.5 * dx)

    k4 = dx * f(y + k3, x + dx)

    return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6.

    if __name__=='__main__':

    t = 0.

    y = 1.

    dt = .1

    ys, ts = [], []

    def func(y, t):

    return t * math.sqrt(y)

    while t <= 10:

    y = runge_kutta(y, t, dt, func)

    t += dt

    ys.append(y)

    ts.append(t)

    exact = [(t ** 2 + 4) ** 2 / 16. for t in ts]

    plt.plot(ts, ys, label='runge_kutta')

    plt.plot(ts, exact, label='exact')

    plt.legend()

    #plt.show()

    结果如下:

    图.3

    示例4:示例1和示例3放在一起进行对比

    import math

    import numpy as np

    import matplotlib.pyplot as plt

    def runge_kutta(y, x, dx, f):

    """ y is the initial value for y

    x is the initial value for x

    dx is the time step in x

    f is derivative of function y(t)

    """

    k1 = dx * f(y, x)

    k2 = dx * f(y + 0.5 * k1, x + 0.5 * dx)

    k3 = dx * f(y + 0.5 * k2, x + 0.5 * dx)

    k4 = dx * f(y + k3, x + dx)

    return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6.

    if __name__=='__main__':

    t = 0.

    y = 1.

    dt = .1

    ys, ts = [], []

    def func(y, t):

    return t * math.sqrt(y)

    while t <= 10:

    y = runge_kutta(y, t, dt, func)

    t += dt

    ys.append(y)

    ts.append(t)

    from scipy.integrate import odeint

    YS=odeint(func,y0=1, t=np.arange(0,10.1,0.1))

    plt.plot(ts, ys, label='runge_kutta')

    plt.plot(ts, YS, label='odeint')

    plt.legend()

    plt.show()

    结果如下

    图.4

    示例5:多个微分方程(欧拉法)

    import numpy as np

    """

    移动方程:

    t时刻的位置P(x,y,z)

    steps:dt的大小

    sets:相关参数

    """

    def move(P, steps, sets):

    x, y, z = P

    sgima, rho, beta = sets

    # 各方向的速度近似

    dx = sgima * (y - x)

    dy = x * (rho - z) - y

    dz = x * y - beta * z

    return [x + dx * steps, y + dy * steps, z + dz * steps]

    # 设置sets参数

    sets = [10., 28., 3.]

    t = np.arange(0, 30, 0.01)

    # 位置1:

    P0 = [0., 1., 0.]

    P = P0

    d = []

    for v in t:

    P = move(P, 0.01, sets)

    d.append(P)

    dnp = np.array(d)

    # 位置2:

    P02 = [0., 1.01, 0.]

    P = P02

    d = []

    for v in t:

    P = move(P, 0.01, sets)

    d.append(P)

    dnp2 = np.array(d)

    """

    画图

    """

    from mpl_toolkits.mplot3d import Axes3D

    import matplotlib.pyplot as plt

    fig = plt.figure()

    ax = Axes3D(fig)

    ax.plot(dnp[:, 0], dnp[:, 1], dnp[:, 2])

    ax.plot(dnp2[:, 0], dnp2[:, 1], dnp2[:, 2])

    plt.show()

    结果如下:

    图.5

    参考:

    展开全文
  • 微分方程求解--补充知识,举例说明如何定义函数文件,高阶常微分方程变量替换 化为一阶常微分方程
  • 常微分方程的数值求解

    千次阅读 2018-04-13 11:59:32
    首先理解一下什么是常微分方程,简单的说就是只有一个未知数的微分方程,具体定义如下: 凡含有参数,未知函数和未知函数导数 (或微分) 的方程,称为微分方程,有时简称为方程,未知函数是一元函数的微分方程称作...
  • Python求解常微分方程

    2021-02-03 16:31:58
    sympy、numpy、scipy、matplotlib是强大的处理数学问题的库,可以执行积分、求解常微分方程、绘图等功能,其开源免费的优势可以与MATLAB媲美。 一阶常微分方程 from sympy import * f = symbols('f', cls=Function...
  • 参考《常微分方程》第三版(王高雄)常微分方程 王高雄 第五章 线性微分方程组_哔哩哔哩 (゜-゜)つロ 干杯~-bilibili​www.bilibili.com前面介绍了线性微分方程的求解,现在考虑线性微分方程组的求解。5.1 存在唯一...
  • 附:微分方程的分类:常微分方程和偏微分方程。1、常微分方程(ODE)是指微分方程的自变量只有一个的方程。最简单的常微分方程,未知数是一个实数...偏微分方程的阶数定义类似常微分方程,但更细分为椭圆型、双曲线...
  • Matlab求常微分方程组的解析解

    千次阅读 2020-03-07 21:53:14
    Matlab求常微分方程的解析解 ...一阶常微分方程求解(无初值) 方程:dy/dt=ay 代码: syms y(t) a; %syms 定义方程用的变量和未知常量(可有可无)。y(t)表示y是关于t的函数,求解t。 eqn=diff(y...
  • 数值分析(5)——求解常微分方程问题的定义Euler折线法梯形法改进Euler法4阶R-K方法 承接自数值分析(4)——数值积分 问题的定义 f(x,y)=dydx,y(a)=ηf(x,y)=\frac{dy}{dx},y(a)=\etaf(x,y)=dxdy​,y(a)=η ...
  • 常微分方程自救手册

    2019-06-09 18:43:00
    常微分方程自救手册 一、微分方程求解 1. 恰当方程 定义 \[ P(x,y)dx + Q(x,y)dy = 0 \] 存在可微函数\(\Phi(x,y)\),使得 \[ d\Phi(x,y)=P(x,y)dx+Q(x,y)dy \] 解法 Step 1:验证是否为恰当方程: \[ \frac{\partial{...
  • ode45是求常微分方程的数值解的首选方法。 matlab提供了求常微分方程数值解的函数。当难以求得微分方程的解析解时,可以求其数值解. matlab中求微分方程数值解的函数有七个:ode45,ode23,ode113,ode15s,ode23s,...
  • 一份较为详细的神经常微分方程学习笔记,共计194页的。其主要涉及问题定义、ODE的基础理论、不同的求解方法、实验以及应用。欢迎感兴趣的朋友下载学习
  • ##定义未知量和待定参数amplitude和detuning amplitude=20 detuning=20 #郎之万方程经典形式 def dfun(t,x,y,q,p): dvec=np.zeros((4)) dvec[0]=-0.5*x-(detuning+0.02*q)*y+amplitude dvec[1]=-0.5*y+(detuning...
  • 在a是常量的情况下,我猜您调用了scipy.integrate.odeint(fun, u0, t, args),其中fun的定义与您的问题相同,u0 = [x0, y0, z0]是初始条件,t是要为ODE求解的时间点序列,args = (a, b, c)是要传递给fun的额外参数。...
  • 目录基本方法求解常微分方程的通解求解常微分方程的初边值问题求解常微分方程求解线性常微分方程组一阶齐次线性微分方程组非齐次线性方程组 基本方法 1.dsolve()函数 求解常微分方程的通解 在求通解问题 % syms y...
  • 常微分方程实验(3.2):解的延拓

    千次阅读 2017-10-22 01:48:04
    嗨~那么之前我们已经学习了(ahhh其实是我已经学习了)如何使用逐步逼近法来求解一个较麻烦但是满足利普希茨条件的初值问题的近似解。但是,现在情况是这样的,还记得之前我们的所求的函数是如何定义的吗?哈当然是...
  • ​ 考虑一个定义在[a,b][a,b][a,b]上的二阶常微分方程边值问题: −u¨(t)+q(t)u(t)=f(t),a<t<bu(a)=0,u(b)=0 -\ddot u(t)+q(t)u(t) = f(t),a<t<b\\ u(a) = 0,u(b) = 0 −u¨(t)+q(t)u(t)=f(t),a<t<...
  • 1.我们通过调用ODE32函数来求解ODE:  [t,y] = ode23('func_name', [start_time, end_time], y(0))  ode45函数使用更高阶的Runge-Kutta公式。  首先我们定义函数,我们创建一个.m文件,输入下面的内容。  ...
  • 自由项f(x)为定义在区间I上的连续函数,即y''+py'+qy=0时,称为二阶系数齐次线性微分方程。若函数y1和y2之比为常数,称y1和y2是线性相关的;若函数y1和y2之比不为常数,称y1和y2是线性无关的。特征方程为:λ^2+p...
  • 练习代码。使用四阶龙格库塔法求解常微分方程,分为了用户定义步长和自动变步长两种。并有一个使用案例,这两种算法的结果与ode45对比。
  • 微分方程-特征线法(转)

    万次阅读 2012-06-03 10:55:57
    数学中的特征线法是求解偏微分方程的一种方法,适用于准...两组常微分方程中的一组用于定义特征线,另一组用以描述解沿给定特征线变化。 设所需求解的准线性偏微分方程为  )(1) 其中 。 取某变量 ,令
  • Python解微分方程

    2018-05-23 15:33:00
    1.求解常微分方程的步骤: from sympy import * init_printing() #定义符号常量x 与 f(x) g(x)。这里的f g还可以用其他字母替换,用于表示函数 x = Symbol('x') f, g = symbols('f g', cls=Function) #用...
  • 讨论了精确初值对刻画方程组解的影响,利用精确初值与关联解之间的关系,定义了模糊微分方程组初值问题解,同时给出了模糊微分方程组的模糊初值问题解存在的判定条件和具体求解方法,以一阶系数模糊线性齐次微分方程组...
  • python对于常微分方程的数值求解是基于一阶方程进行的,高阶微分方程必须化成一阶方程组,通常采用龙格-库塔方法. scipy.integrate模块的odeint模块的odeint函数求常微分方程的数值解,其基本调用格式为: sol=...
  • 注:一阶微分方程求解,在定义函数时,要确保形式为y’=f(t,y),则函数定义时如下例所定义,运算时输入函数即为y’=(y+3*t)/t^2 注:对于二阶及以上则将高阶转为一阶 例: 1.3 刚性问题 例: 利用ode15s函数...
  • Matlab中ode系列函数的用法 Matlab的ode系列函数(统称odesolver)专门... Matlab中采用ode系列函数求解常微分方程(组)的初值问题一般需要两个文件:主程序文件(如:main.m)和定义函数odefun的文件(即odefun.m...
  • 此示例将说明如何定义常微分方程给出的简单模型。 我们将模拟平面数学摆,如图所示。 m是质量,L是从支撑到质心的距离。 让我们假设弦是不可扩展的和无质量的,而且,让我们忽略空气的阻力并假设引力场与g作为...

空空如也

空空如也

1 2 3 4
收藏数 62
精华内容 24
关键字:

常微分方程求解定义