my matlab.rar_instantkme_mumerical algebra_pde

%FDE约束问题2-D： %目标函数 min J(u,y)=\frac{1}{2}iint(y-y1)^{2}dxdt+\frac{\gamma}{2}iint u^{2}dxdt； %约束条件：（^{C}_{0}D_{t}^{\alpha}-^{R}D_{x_1}^{\beta1}-^{R}D_{x_1}^{\beta2})y(x_1,x_2,t)+u(x_1,x_2,t)=f(x_1,x_2,t). %y1:观察项 ; %u:控制项; %% clear; clc; tic; alpha = 0.5; % 时间分数阶 beta1 = 1.5; beta2 = 1.5; % 空间分数阶 gam = 1e-6; % 正则化参数 d = 4; % 网格水平 n = 2^d+2 ; % 空间步数 hx = 1/(n-1); % 空间步长 hy = hx; nt = 2^(d); % 时间步数 ht = hx*hy; % 时间步长 nu = hx*hy*ht; m = n*n*nt; % 期望项 tt = ones(nt,1); %时间节点取one(1,nt),由于观察函数无t. x = hx*(1:1:n)'; %空间节点 [x1,x2,tt1] = ndgrid(x,x,tt); Y = 10*cos(10*x1).*sin(x1.*x2).*tt1; %函数无t,生成三维张量矩阵后拉伸 % u0 = Y(:,:,1); % figure(1) % surf(u0); % Y = permute(Y,[2,1,3]); %变换维度 y0 = reshape(Y,m,1); %观测函数,拉伸成向量 %% 构造 % [M1 0 A'][y]=[f1] % [0 gam*M2 M3'][u]=[f2] % [A M3 0 ][p]=[f3] % 由于积分使用梯形矩形准则，则M1，M2,M3为对角阵，取单位阵的倍数 % M1 = m1*I; M2 = m2*I; M3 = I. % A = kron(Ixy,Dt)-kron(Lspace,It) %=============================================================== m1 = nu; m2 = gam*nu; m3 = nu; f1 = m1*y0; f3 = 0*f1; % f2 = f3; % 分数阶 A 的生成，向量形式存储 % 空间 Toeplitz 列向量 v_beta1 = cumprod([1, 1 - ((beta1+1) ./ (1:n))]); L1c = v_beta1(2:end)'; L1c(1) = L1c(1) + v_beta1(2); L1c(2) = L1c(2) + v_beta1(1); L1c = -0.5*hx^(-beta1)*L1c; %左右导系数相等,生成对称的Toeplitz矩阵 v_beta2 = cumprod([1, 1 - ((beta2+1) ./ (1:n))]); L2c = v_beta2(2:end)'; L2c(1) = L2c(1) + v_beta2(2); L2c(2) = L2c(2) + v_beta2(1); L2c = -0.5*hy^(-beta2)*L2c; % 时间 v_alpha = cumprod([1, 1 - ((alpha+1) ./ (1:nt))]); Cc = ht^(-alpha)*v_alpha(1:nt)'; Cr = ht^(-alpha)*[v_alpha(1);zeros(nt-1,1)]; % 生成矩阵块中A,BTTB矩阵行列信息 t1 = nu*gen3d(L1c,L2c,Cc,Cr,n,n,nt); t2 = nu*gen3d(L1c,L2c,Cr,Cc,n,n,nt); % 循环矩阵的特征值(FFT变换) tev1 = fftn(t1); tev2 = fftn(t2); % Ax = b; % @(x)afun(m1,m2,m3,tev1,tev2,x); % 生成近似的循环矩阵的特征值（FFT） ev1 = genl1ev3d(t1); ev2 = genl1ev3d(t2); %近似矩阵 C Cy=x; %@(x)pfun(ev1,ev2,x) tol = 1e-6; maxit = 1000; b = [f1;f3;f3]; C = toeplitz(Cc,Cr); L1 = toeplitz(L1c); L2 = toeplitz(L2c); L = kron(L1,eye(n))+kron(eye(n),L2); A = kron(C,eye(n*n))-kron(eye(nt),L); A = nu*sparse(A); O = sparse(m,m); I = speye(m); M1 = m1*I; M2 = m2*I; M3 = m3*I; M = [M1 O A';O M2 M3';A M3 O]; y0 = M*b; y1 = afun(m1,m2,m3,tev1,tev2,m,b); % ynum = minres(@(x)afun(m1,m2,m3,tev1,tev2,m,x),b,tol,maxit,... % @(x)pfun(ev1,ev2,m1,m2,gam,m,x)); % ynum = gmres(@(x)afun(m1,m2,m3,tev1,tev2,m,x),... % b,[],tol,maxit);%,@(x)pfun(ev1,ev2,m1,m2,gam,m,x)); % % ynum = pcg(@(x)afun(m1,m2,m3,tev1,tev2,m,x),b,tol,maxit);%,@(x)pfun(ev1,ev2,m1,m2,gam,m,x)); % y = ynum(1:m); % err = norm(y-y0)/norm(y)