分段偏微分方程求解！

机器视觉

下面是我编的一个关于分段偏微分方程的程序，由于初学问题很多，希望各位高手指点一下，问题在那里？急用！
function y=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]内的划分节点数

% 差分计算的初始值
m=round(T/k)+1;
r=c*k/h;
r2=r^2;
r22=r^2/2;
s1=1-r^2;
s2=2-2*r^2;

%情况一
if x>=0&x<=L2 %Compute first and second rows
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
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

U=U';

这个函数的输入初值为：
>> 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;

在此先谢谢了！

happy

x的初始值呢？

机器视觉

我求解的函数中x的值是一个区间：0<=x<=3.5(单位：m)

机器视觉

自己顶一顶！那位高手帮忙看看！

happy

你这个程序专业化有点太强，个别地方看不懂，不过程序中有一些比较低级的错误请先自己运行根据提示修改一下

比如

1. r=c*k/h;

复制代码

这个语句中的c就没有事先赋值

机器视觉

谢谢！