更正：分段偏微分方程求解

机器视觉

我的前一贴里的程序里的c值是个变量，原来的程序有些问题，重新发上来，希望高手们帮我改改，问题在那里？为了便于了解原方程组，在附件里列出了要解的分段偏微分方程组：
function U=finedif(x,M,g,A2,d1,d2,D1,D2,L1,L2,Eg,E,u,vm,P0,p0,T,h,k)

% finedif - 有限差分法求解波动方程的函数
% M - 下落岩体物体总重量,kg
% g - 重力加速度,m/s2
% A2- 活柱体的横截面积,m2
% d1,d2 - 活柱体的内径及外径,m
% D1,D2 - 油缸体的内径及外径,m
% L1 - 活柱体的长度,m
% L2 - 油缸体受高压油作用段的长度,m
% Eg - 液体自身体积弹性模量,MPa
% E - 活柱体及油缸体材料的弹性模量,MPa
% u - 活柱体及油缸体材料的泊松比
% vm - 下落物体的冲击速度,m/s
% P0 - 液体的初始压力,MPa
% p0 - 液体平衡状态的密度,kg/m3

% T - 有限差分计算的时间
% h - 有限差分长度的步长
% k - 有限差分时间的步长

L=L1+L2;

% L - 单体液压支柱的总长

if x>=0&x<=L2
    Ec2=1/((1/Eg)+(2/E)*((D2^2+D1^2)/(D2^2-D1^2)+u));
    c2=sqrt(Ec2/p0);
elseif x>=L2&x<=L
    Ec1=1/((1/Eg)+(2/E)*((d2^2+d1^2)/(d2^2-d1^2)+u));
    c1=sqrt(Ec1/p0);
end

% Ec1 - 考虑活柱体径向变形的液体有效体积弹性模量,MPa
% Ec2 - 考虑油缸体径向变形的液体有效体积弹性模量,MPa
% c1 - 压力波在活柱体内液体中的传播速度,m/s
% c2 - 压力波在油缸体内液体中的传播速度,m/s

if x>=0&x<=L2
    f=(M*g/A2-P0)*x/Ec2;
elseif x>=L2&x<=L
    f=(M*g/A2-P0)*L2/Ec2+(M*g/A2-P0)*x/Ec1;
end

% f - 波动方程初始条件u(x,0)=f(x)

if x>=0&x<L
    g=0;
elseif x>=L
    g=vm;
end

% g - 波动方程初始条件u'(x,t)|t=0=g(x)

if x>=0&x<=L2
    n1=round(L2/h)+1;
elseif x>=L2&x<=L
    n2=round((L-L2)/h)+1;
end

% n1 - 在区域[0,L2]内的划分节点数
% n2 - 在区域[L2,L]内的划分节点数

%情况一
if x>=0&x<=L2 %Compute first and second rows
% 差分计算的初始值
m=round(T/k)+1;
r=c2*k/h;
r2=r^2;
r22=r^2/2;
s1=1-r^2;
s2=2-2*r^2;
for i=2:(n1-1)
    U(i,1)=feval(f,h*(i-1));
    U(i,2)=s1*feval(f,h*(i-1))+k*feval(g,h*(i-1))...
        +r22*(feval(f,h*i)+feval(f,h*(i-2)));
end
for j=3:m %Compute remaining rows
    for i=2:(n1-1)
        U(i,j)=s2*U(i,j-1)+r2*(U(i-1,j-1)+U(i+1,j-1))-U(i,j-2);
    end
end
%情况二
elseif x>=L2&x<=L %Compute remaining rows
% 差分计算的初始值
m=round(T/k)+1;
r=c1*k/h;
r2=r^2;
r22=r^2/2;
s1=1-r^2;
s2=2-2*r^2;
for i=2:(n1-1)
    U(i,1)=feval(f,h*(i-1));
    U(i,2)=s1*feval(f,h*(i-1))+k*feval(g,h*(i-1))...
        +r22*(feval(f,h*i)+feval(f,h*(i-2)));
end
for j=3:m %Compute remaining rows
    for i=2:(n1-1)
        U(i,j)=s2*U(i,j-1)+r2*(U(i-1,j-1)+U(i+1,j-1))-U(i,j-2);
    end
end
end
输入值：
>> M=5000;g=9.8;a2=0.0055;d1=0.084;d2=0.096;D1=0.10;D2=0.114;L1=1.912;L2=1.588;
>> Eg=2100;E=210000;u=0.3;vm=3;P0=30;p0=1000;T=0.010;h=0.1;k=0.01;A2=0.0055;
>> x=0:0.1:3.5;
>> U=finedif(x,M,g,A2,d1,d2,D1,D2,L1,L2,Eg,E,u,vm,P0,p0,T,h,k)
谢谢！！

AaronSpark

错误好多阿

首先你程序里的x是数值，你调用的时候却又输入一组向量

第二个问题feval用法不对，先看看feval的帮助吧

其他的还没看

机器视觉

谢谢楼上的热心指教。
我参照了一下《数值方法Matlab版》中关于差分方法解偏微分方程的方法，程序中可以这样使用变量x;
程序里引用函数feval的部分也是参照其中的。
但是书中的变量x是在一个连续区间，不是分段区间。所以我在解书中方程的时候完全可以，但自己改成分段就有问题了。

机器视觉

顶一顶！哪位好心人帮忙看看

happy

[y1,y2,...] = feval(fhandle,x1,...,xn)
[y1,y2,...] = feval(function,x1,...,xn)

happy

至于x问题，你可以采用循环调用，计算不同x下的结果