标签: differential-equations

我有一组定义3D轮廓的3D点.我想要做的是获得与此轮廓相对应的最小曲面表示(请参阅维基百科中的最小曲面).基本上这需要求解非线性偏微分方程.

在Matlab中,使用该pdenonlin函数几乎是直截了当的(参见Matlab的文档).可以在此处找到用于解决最小表面问题的示例:单位磁盘上的最小表面问题.

我需要在Python中实现这样的实现,但我知道我还没有找到任何关于如何实现这一目标的Web资源.

任何人都可以指出这种实施的任何资源/例子吗？

谢谢,米格尔.

UPDATE

我想要找到的3D表面(理想情况下是三角形网格表示)受这组3D点的限制(如图所示,这些点位于最佳拟合平面中):

好吧,做一些研究我发现这个最小的表面问题与双调和方程的解决方案有关,我还发现薄板样条是这个方程的基本解.

所以我认为这种方法是尝试使用薄板样条来拟合表面的稀疏表示(由点的3D轮廓给出).我在scipy.interpolate中找到了这个例子,其中使用薄板样条插值散射数据(x,y,z格式)以获得均匀网格(XI,YI)上的ZI坐标.

出现两个问题:(1)薄板样条插值是否是从三维轮廓点集计算表面问题的正确方法？(2)如果是这样,如何使用NON-UNIFORM网格在scipy上执行薄板插值？

再次感谢!米格尔

更新:MATLAB中的实现(但它不适用于SCIPY PYTHON)

我使用Matlab的函数跟踪了这个例子,tpaps并在均匀网格上获得了适合我轮廓的最小曲面.这是Matlab的结果(看起来很棒!):

但是我需要在Python中实现它,所以我使用包scipy.interpolate.Rbf和thin-plate函数.这是python中的代码(XYZ包含轮廓中每个点的3D坐标):

GRID_POINTS = 25
x_min = XYZ[:,0].min()
x_max = XYZ[:,0].max()
y_min = XYZ[:,1].min()
y_max = XYZ[:,1].max()
xi = np.linspace(x_min, x_max, GRID_POINTS)
yi = np.linspace(y_min, y_max, GRID_POINTS)
XI, YI = np.meshgrid(xi, yi)

from scipy.interpolate import Rbf
rbf = Rbf(XYZ[:,0],XYZ[:,1],XYZ[:,2],function='thin-plate',smooth=0.0)
ZI = rbf(XI,YI)

然而,这是结果(与Matlab中获得的结果完全不同):

很明显,scipy的结果与最小的表面不对应. …

我将动态模型设置为ODE的(僵硬)系统.我目前用CVODE解决了这个问题(来自Assimulo python包中的SUNDIALS包),一切都很好.

我现在想要为问题添加一个新的3D散热器(具有与温度相关的热参数).我的想法是使用现有的FEM或FVM框架为我提供一个界面,使我可以轻松地为(t, y)常规提供3D块,并获得,而不是从头开始为3D热方程编写所有方程式.支持残差y'.原理是使用FEM系统中的方程而不是求解器.CVODE可以利用稀疏性,但预计组合系统的解决速度将慢于FEM系统自身解决的速度,并为此量身定制.

# pseudocode of a residuals function for CVODE
def residual(t, y):

    # ODE system of n equations 
    res[0] = <function of t,y>;
    res[1] = <function of t,y>;
    ...
    res[n] = <function of t,y>;

    # Here we add the FEM/FVM residuals
    for i in range(FEMcount):
        res[n+1+i] = FEMequations[FEMcount](t,y)

    return res

我的问题是(a)这种方法是否合理,(b)是否有一个FEM或FVM库可以轻松地让我将其视为一个方程组,这样我就可以"将它"粘贴到我现有的一套ODE方程.

如果不能让两个系统共享同一个时间轴,那么我将不得不以步进模式运行它们,我在那里运行一个模型一小段时间,更新另一个模型的边界条件,运行那个,更新第一个模型的BC,等等.

我对精彩的图书馆FiPy有一些经验,我期望最终以上述方式使用该库.但我想知道其他系统在这种性质问题上的经验,以及我错过的其他方法.

编辑:我现在有一些看似有效的示例代码,显示了如何使用CVODE解决FiPy网格扩散残差.然而,这只是一种方法(使用FiPy),其余的问题和疑虑仍然存在.欢迎任何建议.

from fipy import *
from fipy.solvers.scipy import DefaultSolver
solverFIPY = DefaultSolver()

from assimulo.solvers import CVode as solverASSIMULO
from …

我正在使用DE系统,我想知道哪个是最常用的python库来解决微分方程(如果有的话).

我的方程是非线性一阶方程.

更具体地说,我对8阶Dormand-Prince嵌入式方法感兴趣,它基于Runge-Kutta和刚性方程.

我使用Numerical Recipes 3但我经常在编译它们的库时遇到问题.我想知道其他选择.

我想在另一个函数中使用函数的导数.如何在Maxima中完成？

例如:

f(x) := 2*x^4;
g(x) := diff(f(x),x)-8;

现在按预期g(x)收益8x^3-8,但是 g(0)给出错误,因为diff(f(0),0)没有意义.但那我该如何正确定义g呢？

我如何在Python中用数字解决ODE？

考虑

\ddot{u}(\phi) = -u + \sqrt{u}

具备以下条件

u(0) = 1.49907

和

\dot{u}(0) = 0

有约束

0 <= \phi <= 7\pi.

最后,我想生成一个参数图,其中x和y坐标是作为u的函数生成的.

问题是,我需要运行odeint两次,因为这是二阶微分方程.我试过让它在第一次运行后再次运行,但它回来时出现雅可比错误.必须有一种方法可以同时运行它两次.

这是错误:

odepack.error:该函数及其雅可比行列式必须是可调用函数

以下代码生成的代码.有问题的行是sol = odeint.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from numpy import linspace


def f(u, t):
    return -u + np.sqrt(u)


times = linspace(0.0001, 7 * np.pi, 1000)
y0 = 1.49907
yprime0 = 0
yvals = odeint(f, yprime0, times)

sol = odeint(yvals, y0, times)

x = 1 / sol * …

我正在尝试开发物理模拟,我想实现一个四阶辛积分方法.问题是我必须弄错数学,因为在使用辛积分器时我的模拟根本不起作用(与模拟工作相当好的四阶Runge-Kutta积分器相比).我一直在谷歌搜索这个,我能找到的只是关于这个主题的科学文章.我试图改编文章中使用的方法,但我没有运气.我想知道是否有人有使用辛积分器的模拟的源代码,最好是模拟引力场,但任何辛积分器都可以.源代码的语言并不重要,但我会欣赏使用C风格语法的语言.谢谢!

我试图在java中创建一个简单的SIR-epidemics模型模拟程序.

基本上,SIR由三个微分方程组定义:
S'(t)= - l(t)*S(t)
I'(t)= l(t)*S(t) - g(t)*I(t)
R'(t)= g(t)*I(t)

S - 易感人群,我 - 感染了人,R - 恢复了人.

l(t)= [c*x*I(t)]/N(T)

c - 接触次数,x - 感染性(与病人接触后生病的概率),N(t) - 总人口(不变).

如何在Java中解决这些微分方程？我不认为我知道有任何有用的方法,所以我的实现会产生垃圾.

public class Main {
public static void main(String[] args) {
    int tppl = 100;
    double sppl = 1;
    double hppl = 99;
    double rppl = 0;
    int numContacts = 50;
    double infectiveness = 0.5;
    double lamda = 0;
    double duration = 0.5;
    double gamma = 1 / duration;
    for (int i = …

我真的不知道从哪里开始这个问题,因为我没有太多经验,但需要使用计算机解决项目的这一部分.

我有一个二阶ODE,它是:

m = 1220

k = 35600

g = 17.5

a = 450000

和b在1000和10000之间,增量为500.

x(0)= 0 

x'(0)= 5


m*x''(t) + b*x'(t) + k*x(t)+a*(x(t))^3 = -m*g

我需要找到最小的b,这样解决方案永远不会是积极的.我知道图形应该是什么样的,但我只是不知道如何使用odeint来获得微分方程的解.这是我到目前为止的代码:

from    numpy    import    *    
from    matplotlib.pylab    import    *    
from    scipy.integrate    import    odeint

m = 1220.0
k = 35600.0
g  = 17.5
a = 450000.0
x0= [0.0,5.0]

b = 1000

tmax = 10
dt = 0.01

def fun(x, t):
    return (b*x[1]-k*x[0]-a*(x[0]**3)-m*g)*(1.0/m)
t_rk = arange(0,tmax,dt)   
sol = odeint(fun, x0, t_rk)
plot(t_rk,sol)
show()

哪个并没有真正产生任何东西.

有什么想法吗？谢谢

如果有人可以帮助解决以下问题,我将不胜感激.我有以下ODE:

dr/dt = 4*exp(0.8*t) - 0.5*r   ,r(0)=2, t[0,1]       (1)

我用两种不同的方式解决了(1).通过Runge-Kutta方法(第4顺序)和ode45Matlab中的方法.我将这两个结果与分析解决方案进行了比较,分析解决方案由下式给出:

r(t) = 4/1.3 (exp(0.8*t) - exp(-0.5*t)) + 2*exp(-0.5*t)

当我根据确切的解决方案绘制每个方法的绝对误差时,我得到以下结果:

对于RK方法,我的代码是:

h=1/50;                                            
x = 0:h:1;                                        
y = zeros(1,length(x)); 
y(1) = 2;    
F_xy = @(t,r) 4.*exp(0.8*t) - 0.5*r;                   
for i=1:(length(x)-1)                              
    k_1 = F_xy(x(i),y(i));
    k_2 = F_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);
    k_3 = F_xy((x(i)+0.5*h),(y(i)+0.5*h*k_2));
    k_4 = F_xy((x(i)+h),(y(i)+k_3*h));
    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
end

并为ode45:

tspan = 0:1/50:1;
x0 = 2;
f = @(t,r) 4.*exp(0.8*t) - 0.5*r;
[tid, y_ode45] …