| 查看: 616 | 回复: 2 | |||
[交流]
求助，matlab 四阶龙格库塔结合二乘法的问题 已有1人参与
|
|
function k format long clear all clc x0 = 0; k0 = 0; lb = 0.1; ub = 0.5; tspan=linspace(0,400,800); [m,n]=ode45(@example,tspan,0); tspan = m; data=n; plot(tspan,data); yexp = data'; options = optimset('Display','iter','TolFun',1e-10,'TolX',1e-10) [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,k0,lb,ub,options,tspan,x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk = %.9f ± %.9f\n',k,ci(1,2)-k) function f = ObjFunc(k,tspan,x0,yexp) % 目标函数 [tspan,Xsim] = ode45(@example2,tspan,0,[],k); Ysim=Xsim'; Ysim=Ysim(2:end); f = Ysim(:,1)-yexp(:,1); function dy = example(t,y) pi=3.14; dy=0.1*y-0.1*y.^3+0.3*sin(2*pi*0.01*t)+sqrt(2)*randn(size(t)); function dy = example2(t,y,k) %ode 模型 %UNTITLED4 Summary of this function goes here % Detailed explanation goes here pi=3.14; dy=0.1*y-0.1*y.^3+k*sin(2*pi*0.01*t)+sqrt(2)*randn(size(t)); 这是一个正确的程序。但是我想问，如果把sqrt(2)*randn(size(t))改为一个已知的1*800数组，然后对应时间t依次代进去，怎么弄？ |
回复此楼» 猜你喜欢
浙大物院表面物理课题组长期招收普博及直博学生
已经有259人回复
-
【征稿】Springer Nature热电专题 Discover Materials / Discover Chemistry
已经有47人回复
-
物理学I论文润色/翻译怎么收费?
已经有135人回复
-
25届硕士申请，材料计算方向
已经有2人回复
-
课题组招1-2名硕士
已经有2人回复
-
机时越来越贵之后，我把一部分 DFT 流程搬到了 AI 计算平台上
已经有0人回复
-
斯坦福大学与Elsevier世界前2%学者新排名查询
已经有1人回复
-
27申博推荐：南京航空航天大学国际前沿院光学方向招收博士生！
已经有19人回复
-
华南师范大学+双一流学科物理学+广东省杰青团队2027年博士招生+机器学习驱动材料设计
已经有1人回复
-
apl投稿状态
已经有2人回复
-
华南师范大学+双一流学科物理学+广东省杰青团队2027年博士招生+机器学习驱动材料设计
已经有1人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 求助用matlab编一个程序，采用龙格库塔发求求解动力学方程的参数问题
已经有5人回复
- 急！求大神用Matlab四阶龙格库塔解个方程！
已经有5人回复
- 请问什么是分段龙格库塔的方法？需要用matlab软件采用这种方法求二阶微分方程的解
已经有3人回复
- Matlab用龙格库塔法求二阶非线性方程
已经有13人回复
- 求matlab 用四阶龙格-库塔法求解微分方程
已经有4人回复
- 龙格库塔法 matlab
已经有6人回复
- 【求助】Matlab中利用四阶龙格-库塔法求解微分方程！！！！
已经有9人回复
月只蓝
主管区长 (职业作家)
-
专家经验: +1059 - 计算强帖: 8
- 应助: 1712 (讲师)
- 贵宾: 8.888
- 金币: 68182.7
- 散金: 1938
- 红花: 443
- 沙发: 4
- 帖子: 4373
- 在线: 3292.4小时
- 虫号: 1122189
- 注册: 2010-10-14
- 专业: 宇宙学
- 管辖: 计算模拟区
2楼2014-10-28 08:37:59
|
function k format long clear all clc x0 = 0; k0 = 0; lb = 0.1; ub = 0.5; tspan=linspace(0,400,800); [m,n]=ode45(@example,tspan,0); tspan = m; data=n; plot(tspan,data); yexp = data'; options = optimset('Display','iter','TolFun',1e-10,'TolX',1e-10) [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,k0,lb,ub,options,tspan,x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk = %.9f ± %.9f\n',k,ci(1,2)-k) function f = ObjFunc(k,tspan,x0,yexp) % 目标函数 [tspan,Xsim] = example2(tspan,k,x0); Ysim=Xsim'; Ysim=Ysim(2:end); f = Ysim(:,1)-yexp(:,1); function dy = example(t,y,y0) pi=3.14; dy=0.1*y-0.1*y.^3+0.3*sin(2*pi*0.01*t)+sqrt(2)*randn(size(t)); function [t,y] = example2(t,k,y0) %ode 模型 fs=4;%采样频率 %T=1/fs;%采样时间 N=4096;%数据点数 %t=(0:N-1)*T;%时间序列 pi=3.14; t_1=linspace(0,400,800); [m1,n1]=ode45(@example,t_1,0); subplot(3,2,1) plot(m1,n1); NFFT=2^nextpow2(N);%求得最接近总采样点的2^n, Y2=fft(n1,NFFT)/N; y2=2*abs(Y2(1:NFFT/20+1)); f2=fs/2*linspace(0,0.1,NFFT/20+1);%频率序列 subplot(3,2,2) plot(f2,y2); %第二个经过SR的频谱图 LPF = ones(4096,1); LPF(18:24)=zeros(7,1); %截止频率100Hz %LPF(14:18)=zeros(5,1); %截止频率100Hz Ylv4 = Y2.*LPF; f4=fs/2*linspace(0,0.1,NFFT/20+1); ylv4=2*abs(Ylv4(1:NFFT/20+1)); subplot(3,2,4) plot(f4,ylv4); %第四个图是幅值置0后对比图 T5=1/N; %t5=(0:N-1)*T5;%时间序列 t5=linspace(0,400,800); ylv5 = ifft(Ylv4*N); ylvlv=ylv5(1:5:3996); ylvlv2=ylvlv'; subplot(3,2,5) plot(4*t5 ,ylv5(1:5:3996)); %画出第个五图 %上面的这些为了得到一个1*800的数据ylvlv2,已经确认得到就是1*800的矩阵，但是用for以后就不能调用ode45，结果出来就是k未定义 len=length(ylvlv2); ty=linspace(0,400,len); y=zeros(1,len); y(1)=a*y0-b*y0^3+ylvlv2(1)+k*sin(2*pi*0.01*ty(1)); for i=1:len-1 h=ty(i+1)-ty(i); q1=h*(a*y(i)-b*y(i)^3+ylvlv2(i+1)+k*sin(2*pi*0.01*ty(i+1))); q2=h*(a*(y(i)+h*k1/2)-b*(y(i)+h*k1/2)^3+ylvlv2(i+1)+k*sin(2*pi*0.01*(ty(i+1)+h/2))); q3=h*(a*(y(i)+h*k2/2)-b*(y(i)+h*k2/2)^3+ylvlv2(i+1)+k*sin(2*pi*0.01*(ty(i+1)+h/2))); q4=h*(a*(y(i)+h*k3/2)-b*(y(i)+h*k3)^3+ylvlv2(i+1)+k*sin(2*pi*0.01*(ty(i+1)+h))); y(i+1)=y(i)+(q1+2*q2+2*q3+q4)/6; end 能再帮忙看看出现什么问题么 灰常感谢 急着用程序 |
3楼2014-10-28 10:06:47
