㈠ matlab 龙格库塔解的复杂微分方程组,自己编的计算简单的方程可以,解复杂的方程组可能不对。请高手指点
用matlab的龙格-库塔法求解复杂微分方程组的关键,在于建立自定义微分方程组函数。你的问题主要是没有利用数组格式来写,没有理顺其之间的关系。具体应按下列格式书写
function dy=odefun(t,y)
ft=0.005*t^3-350*t^2+0.002*t;
hy1=0.5*y(1)^3+0.002*y(1)^2+0.005;
dy(1)=500*ft*sin(t)/(hy1*y(1));
dy(2)=y(2)*(y(1)+(700-2*y(2))*tan(0.5*t))/(700-y(1)*tan(0.5*t));
dy=dy(:);
end
调用格式为
y0=[0.001,0];
[t,y]=ode45(@odefun,[0.58:-0.001:0],y0)
运行结果
㈡ matlab龙格库塔法解微分方程
function [Y] = RK45(t,X,f,h)
K1=f(t,X);
K2=f(t+h/2,X+h/2*K1);
K3=f(t+h/2,X+h/2*K2);
K4=f(t+h,X+h*K3);
Y=X+h/6*(K1+2*K2+2*K3+K4);
end
以上是4阶龙格库塔法的代码:
自己写函数,存为f.m
function dxdt = f (t,x)
dxdt(1)=exp(x(1)*sin(t))+x(2);
dxdt(2)=exp(x(2)*cos(t))+x(1); % x(1)是你的f,x(2)是你的g
dxdt=dxdt(:);
end
自己给出t0,x0,h的值(初始时间,初值,步长)
如果求t0到t1的轨迹的话:给个例子如下
t0=0;t1=5;h=0.02;x0=[-1;-1];
T=t0:h:t1;X=zeros(length(x0),length(T));X(:,1)=x0;
for j=1:length(T)-1
X(:,j+1)=RK45(T(j),X(:,j),@(t,x) f(t,x),h);
end
plot(T,X(1,:));
hold on;
plot(T,X(2,:),'r');
具体参数自己设置
㈢ 用 四阶龙格库塔 解三阶微分方程 公式和程序
初值给的不够啊。
不仅要给y,还要给y的一阶导和二阶导。否则数值解没法弄。
理论解的话,可能还可以含有一些系数。
==============================================================
在Matlab下输入:edit,然后将下面两行百分号之间的内容,复制进去,保存
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function DYDt=_amanm(t,Y)
y=Y(1);
dydt=Y(2);
d2ydt2=Y(3);
DYDt=[Y(2);Y(3);(10-10*y-dydt-0.11d2ydt2)/0.001];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
在Matlab命令行下面输入:
t_start=0;
t_end=1;
yy0=[0.1;0;0]; %初值
[t,y]=ode45('_amanm',[0,t_end],yy0);
plot(t,y(:,1));
xlabel('t');
ylabel('y')
y(end,1)
得到的结果:
ans =
1.0056
用昨天给的C代码计算出来的结果是:
1.005630
㈣ mathematic用龙格库塔法解二阶方程
……所谓龙格库塔法,通俗地说,就是把一个n阶的常微分方程,整理成n个形如 f'(t)=g(t,f(t)) (注意此时右侧不含 f(t) 的导数)的一阶常微分方程组再加以求解的方法。你的方程整理成龙格库塔所需要的形式就是:
x'[t] = y[t]
y'[t] = (-c y[t] - k x[t] - F[t])/m
再考虑到你找到的这段代码本身对格式的要求,就有:
m=1;c=1;k=1;f=1;
stepsize=0.01;t=1;
sol=NestList[
RungeKutta[{y,(-cy-kx-f)/m},#,{x,y},stepsize]&,{0.,0.},Round[t/stepsize]];
剩下的你自己弄吧……算了,考虑到你的方程里的F[t]涉及了t的具体值,又多了一点点坎(在要使用你找到的这个函数的前提下),再来一个F[t_]=Sin[t]的代码示例吧:
Clear[f,t]
m=1;c=1;k=1;
f[t_]=Sin[t];
stepsize=0.01;tmin=0;tmax=1;
sol=FoldList[
RungeKutta[{y,(-cy-kx-f[#2])/m},#,{x,y},
stepsize]&,{0.,0.},Range[tmin,tmax,stepsize]];
㈤ 龙格库塔的文字介绍
龙格-库塔(Runge-Kutta)法
到目前为止,我们已经学习了多步法,例如:亚当斯-巴什福思(Adams
-Bashorth)法,亚当斯-莫尔顿(Adams-Monlton)法,都是常微分方程的积分
方法。它们需要在每一次迭代时重新计算一遍等式右边的结果(非线性隐含问
题忽略计算多个 f (ω)值的可能性)
龙格-库塔(Runge-Kutta)法是一种不同的处理,作为多级方法为人们所知。
它要求对于一个简单的校正计算多个 f 的值。
下面,我们列出了 3 种最流行的龙格-库塔(Runge-Kutta)法:
改进的欧拉方法(精度:p=2):
V a = V n + Δtf (V n,tn)
2
Δt)二阶格式
V n+1 = V n +Δtf (V a,tn +
2
Hevn’s 方法(p=2):
这是另一种二阶格式:
V a = V n +Δtf (V n,tn)
V n = V n +
+1 Δt[ f (V n,tn) + f (V a,tn +Δt)]
2
注意: f (Vn,tn)在运算中应该只被计算一次。
四次龙格-库塔(Runge-Kutta)法(p=4):
这是一个 4 阶格式。这次我们写的形式有点不同:
a = Δtf (V n,tn)
b = Δtf (V n + 1 a,tn + 1
2 2 Δt)
c = Δtf (V n + 1 b,tn + Δt)
1
2 2
d = Δtf (V n + c,tn +Δt)
V n =V n +
+1 1 (a +2b +2c +d)。
6
㈥ matlab龙格库塔法求解微分方程,怎么编程
用matlab龙格库塔法(ode45)求解微分方程组的方法。
第一步:建立微分方程阻的自定义函数myodefun(t,z),其主要内容
dz(1)=1.2*z(1)-0.3*z(1)*z(2);
dz(2)=0.05*z(1)*z(2)-0.1*z(2);
第二步:用ode45函数求解。
[t,z]=ode45(@myodefun,[0,30],[10,3])
第三步:用plot函数绘制t-z1(t)和t-z2(t)的图形。
计算结果
㈦ 龙格库塔求解二阶微分方程组的MATLAB编程
MATLAB求解x''+0.7x'+0.8x'|x'|+25.6x-25.6x³=0二阶微分方程组的方法,可以按下列步骤进行:
1、建立自定义函数func()
function
f
=
func(t,x)
%x''+0.7x'+0.8x'|x'|+25.6x-25.6x³=0
f(1)=x(2);
f(2)=25.6*x(1)^3-25.6*x(1)-0.8*x(2)*abs(x(2))-0.7*x(2);
f=f(:);
2、建立龙格库塔算法函数runge_kutta()
调用格式:[t,x]
=
runge_kutta(@(t,x)func(t,x),x0,h,a,b);
3、然后根据x和x'数据,绘制出x(t)、x′(t)的图形。
plot(x(:,1),x(:,2))
㈧ MATLAB中已知系统微分方程及初始值用欧拉法和龙格库塔法解一阶微分方程
function Euler
%欧拉法和龙格库塔算法解一阶常微分方程源代码
%例子dy/dx=-y+x+1
f=inline('-y+x+1','x','y'); %微分方程的右边项
dx=0.5; %x方向步长
xleft=0; %区域的左边界
xright=10; %区域的右边界
xx=xleft:dx:xright; %一系列离散的点
n=length(xx); %点的个数
y0=1;
%%(1)欧拉法
Euler=y0;
for i=2:n
Euler(i)=Euler(i-1)+dx*f(xx(i-1),Euler(i-1));
end
%%(2)龙格库塔法
RK=y0;
for i=2:n
k1=f(xx(i-1),RK(i-1));
k2=f(xx(i-1)+dx/2,RK(i-1)+k1*dx/2);
k3=f(xx(i-1)+dx/2,RK(i-1)+k2*dx/2);
k4=f(xx(i-1)+dx,RK(i-1)+k3*dx);
RK(i)=RK(i-1)+dx*(k1+2*k2+2*k3+k4)/6;
end
%%Euler和Rk法结果比较
plot(xx,Euler,xx,RK)
hold on
%精确解用作图
syms x
rightsolve=dsolve('Dy=-y+x+1','y(0)=1','x');%求出解析解
rightdata=subs(rightsolve,xx);%将xx代入解析解,得到解析解对应的数值
plot(xx,rightdata,'r*')
legend('Euler','Runge-Kutta','analytic')
㈨ 龙格库塔求解微分方程组
只要理解了龙格库塔,这就很容易了。
定义函数,
f=func(x,y)
if y(1)>126,
f=[f1(y(1),y(2)), f2(y(1),y(2))];
else
f=[g1(y(1),y(2)), g2(y(1),y(2))];
end
然后就简单了
㈩ 龙格库塔法的基本原理
该算法是构建在数学支持的基础之上的。对于一阶精度的拉格朗日中值定理有:
对于微分方程:y'=f(x,y)
y(i+1)=y(i)+h*K1
K1=f(xi,yi)
当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值,那么就会得到二阶精度的改进拉格朗日中值定理:
y(i+1)=y(i)+[h*( K1+ K2)/2]
K1=f(xi,yi)
K2=f(x(i)+h,y(i)+h*K1)
依次类推,如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km,并用他们的加权平均数作为平均斜率K*的近似值,显然能构造出具有很高精度的高阶计算公式。经数学推导、求解,可以得出四阶龙格-库塔公式,也就是在工程中应用广泛的经典龙格-库塔算法:
y(i+1)=y(i)+h*( K1+ 2*K2 +2*K3+ K4)/6
K1=f(x(i),y(i))
K2=f(x(i)+h/2,y(i)+h*K1/2)
K3=f(x(i)+h/2,y(i)+h*K2/2)
K4=f(x(i)+h,y(i)+h*K3)
通常所说的龙格-库塔法是指四阶而言的,我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式