求解扩散方程

saintpmj

求解扩散方程

saintpmj

不明白，方程

风花雪月

saintpmj 发表于 2010-10-9 10:03

                               
登录/注册后可看大图

不明白，方程

什么不明白？
没看懂方程还是别的什么意思？
请把你的问题描述清楚

saintpmj

哈哈，我现在已经解出来了，不过，我现在又有新麻烦了

saintpmj

回复 风花雪月 的帖子

%计算每个三角单元
L = cell(600,1);
%b = zeros(600,1);
%c = zeros(600,1);
a = double(zeros(331,331));
for i = 1:1:600
    m = B(i,1);
    n = B(i,2);
    l = B(i,3);
    b(m) = (A(B(i,2),2)-A(B(i,3),2))*3/pi;
    b(n) = (A(B(i,3),2)-A(B(i,1),2))*3/pi;
    b(l) = (A(B(i,1),2)-A(B(i,2),2))*3/pi;
    c(m) = (A(B(i,1),1)-A(B(i,3),1))*3/pi;
    c(n) = (A(B(i,1),1)-A(B(i,3),1))*3/pi;
    c(l) = (A(B(i,2),1)-A(B(i,1),1))*3/pi;
    L(m,1) = {(A(B(i,2),1)*A(B(i,3),2)-A(B(i,3),1)*A(B(i,2),2)+(A(B(i,2),2)-A(B(i,3),2))*x+(A(B(i,1),1)-A(B(i,3),1))*y)*3/pi};
    L(n,1) = {(A(B(i,3),1)*A(B(i,2),2)-A(B(i,1),1)*A(B(i,3),2)+(A(B(i,3),2))-A(B(i,1),2)*x+(A(B(i,1),1)-A(B(i,3),1))*y)*3/pi};
    L(l,1) = {(A(B(i,1),1)*A(B(i,2),2)-A(B(i,2),1)*A(B(i,1),2)+(A(B(i,1),2))-A(B(i,2),2)*x+(A(B(i,2),1)-A(B(i,1),1))*y)*3/pi};
   %三角形计算9次
    for j = 1:1:3
        M = B(i,j);
        for k = 1:1:3
            N = B(i,k);
            %二重积分计算
            syms x y;
            I = int(L{M,1}.*L{N,1},y,0,1-x);
            P = int(I,x,0,1);
            P = eval(P);%类型转换
            %计算刚度矩阵
           a(M,N) = a(M,N)+((b(M).*b(N)+c(M).*c(N)));
        end
    end
end

saintpmj

a(M,N) = a(M,N)+((b(M).*b(N)+c(M).*c(N)));
这一行总是说，维数不匹配

Rainyboy

回复 saintpmj 的帖子

b和c都是一维数组，不适合在已经有了索引M和N的情况下用.*吧。
或者你的意思是（？）：
a(M,N) = a(M,N)+((b(M)*b(N)+c(M)*c(N)));

风花雪月

saintpmj 发表于 2010-10-13 17:27

                               
登录/注册后可看大图

a(M,N) = a(M,N)+((b(M).*b(N)+c(M).*c(N)));
这一行总是说，维数不匹配

代码不全，比如缺少B，无法运行
单纯从代码上看没有看出这里有什么问题

Rainyboy

回复 风花雪月 的帖子

%b = zeros(600,1);
%c = zeros(600,1);
这个注释是不是在说明B？

风花雪月

Rainyboy 发表于 2010-10-15 17:22

                               
登录/注册后可看大图

回复 风花雪月 的帖子

%b = zeros(600,1);

不是的，B肯定是另有来源的
b、c应该是坐标，而B应该是节点编号列表

smtmobly

你这个方程是你自己写成积分形式的还是得到的时候就是这个形式，如果有pde形式，可以写出来嘛？

saintpmj

回复 smtmobly 的帖子

我是在arridge他的论文上看到的，我要解这个方程。

saintpmj

我当时，是把P当成一个数值加到a（m,n）里的，其实，后来，我才知道，P 是一个一维数组，不能当成数值直接与数相加

saintpmj

有谁做过关于光学层析成像方面的论文呀。或者实验，我求解出的总有负值，有负值是不正确的。