python - 将代码从matlab转换为python,出现问题
问题描述
我在这里找到了用matlab编写的反向欧拉的实现, 计算步长的方程是:
function [t,y]=beul(f,df,y0,t0,tf,n)
h=(tf-t0)/n;
t=linspace(t0,tf,n+1);
y=zeros(n+1,length(y0));
y(1,:)=y0;
for i=1:n
x0=y(i,:)’;
x1=x0-inv(eye(length(y0))-h*feval(df,t(i),x0))*(x0-h*feval(f,t(i),x0)’-y(i,:)’);
while (norm(x1-x0)>0.0001)
x0=x1;
x1=x0-inv(eye(length(y0))-h*feval(df,t(i),x0))*(x0-h*feval(f,t(i),x0)’-y(i,:)’);
end
y(i+1,:)=x1’;
end
end
那么这个函数被称为定义函数和雅可比系统:
function yp=volt(t,y)
a=4;
c=1;
yp(1)=a*(y(1)-y(1)*y(2));
yp(2)=-c*(y(2)-y(1)*y(2));
end
function y=dvol(t,x)
a=4;
c=1;
y(1,1)=a*(1-x(2));
y(1,2)=-a*x(1);
y(2,1)=c*x(2);
y(2,2)=-c*(1-x(1));
并且反向欧拉被称为:
[t,y]=beul(’volt’,’dvol’,[2,1],0,10,1000);
我已经翻译了python中的代码:
class Backward(Euler):
def solve(self):
for i in range(len(self.time)-1):
u0 = self.u[i].T
u1 = u0 - np.linalg.inv(np.eye(len(self.dydt.u0)) - self.dt * self.dydt.df(self.time[i+1],u0)) * (u0 - self.dt * self.dydt.f(self.time[i+1], u0).T - self.u[i].T )
error = np.array([1.0])
iters = 0
while True:
try:
u0 = u1.T
u1 = u0 - np.linalg.inv(np.eye(len(self.dydt.u0)) - self.dt * self.dydt.df(self.time[i+1],u0)) * (u0 - self.dt * self.dydt.f(self.time[i+1], u0).T - self.u[i].T )
iters += 1
error = np.abs(u1-u0)
if np.sum(np.abs(error)) <= toll:
break
except ValueError as ex:
print('Error occurred in Implicit-NR Euler Solvers: %s' %ex.args)
return None , None
self.u[i+1] = u1.T
然后我定义了雅可比矩阵如下:
def f(self,ti,ui):
return np.array([function(ti,ui) for function in self.func])
def df(self, t, u, **params):
eps = 1e-12
J = np.zeros([len(u), len(u)], dtype = np.float)
for i in range(len(u)):
u1 = u.copy()
u2 = u.copy()
u1[i] += eps
u2[i] -= eps
f1 = self.f(t, u1, **params)
f2 = self.f(t, u2, **params)
J[ : , i] = (f1 - f2) / (2 * eps)
return J
如果我尝试运行单方程问题,这些方法效果很好(我与其他求解器进行了比较)
但问题是matlab产品的行为不同!所以我不知道如何将产品修复为在 python 中相同,因为当我为系统运行代码时(例如 matlab 解决的相同问题)
eq1 = lambda t,u : a*(u[0]-u[0]*u[1]);
eq2 = lambda t,u : -c*(u[1]-u[0]*u[1]);
func1 = np.array([eq1,eq2])
y0 = np.array([2.,1.])
我收到了这个错误:
Running Newton-Rapson Backward Euler ....
Error occurred in Implicit-NR Euler Solvers: could not broadcast input array from shape (2,2) into shape (2)
那么如何定义与matlab计算相同的产品(它的df也是2x2矩阵)以修复python的方法?
解决方案
我使用 numpy.dot 产品解决了:
u1 = u0 - np.linalg.inv(np.eye(len(self.dydt.u0)) - self.dt * self.dydt.df(self.time[i+1],u0)).dot(u0 - self.dt * self.dydt.f(self.time[i+1], u0).T - self.u[i].T )
推荐阅读
- python - 如何将在多行中具有相同键的字典转换为数据框
- laravel - 试图在多态模型中调用imageable_type的关系
- r - 如何在 R 中根据因子变量的每个值的不同比例,从数据集中抽取与大小成比例的随机样本
- python - cv2.inrange() 用于多种颜色检测(要检测63种不同颜色的情况)
- firebase - Flutter:如何获取多个下载网址?
- dialogflow-es - 每当我尝试使用“立即试用”功能时,Dialogflow 都会自动让我退出
- php - Fullcalendar 未在 Laravel 中呈现
- java - 如何在屏幕关闭背景可播放选项的情况下播放 android 收音机
- bash - 尝试使用选择命令 Bash 通过用户选择来执行功能
- sql - 有没有办法在SQL中找到超过n天的最新日期?